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Abstract 

We have completed the evaluation of all mass-dependent QED contributions to the muon 
g — 2, or a^, in two or more different formulations. Their numerical values have been greatly 
improved by an extensive computer calculation. The new value of the dominant term 
A2 {m^/nie) is 132.6823 (72), which supersedes the old value 127.50 (41). The new value of 
the three-mass term (m^/rn-e, m^/rriT-) is 0.0376 (1). The term A2 {m^/nir) is crudely es- 
timated to be about 0.005 and may be ignored for now. The total QED contribution to is 
116 584 719.58 (0.02)(1.15)(0.85) x 10"i\ where 0.02 and 1.15 are uncertainties in the and 
terms and 0.85 is from the uncertainty in a measured by atom interferometry. This raises the 
Standard Model prediction by 13.9 x 10~^^, or about 1/5 of the measurement uncertainty of a^. It 
is within the noise of current uncertainty (~ 100 x 10^^^) in the estimated hadronic contributions 
to a^. 
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I. INTRODUCTION AND SUMMARY 

The latest measured value of the anomalous magnetic moment of negative muon is Q| 

a^_(exp) = 11 659 214 (8) (3) x 10"^° (0.7 ppm), (1) 

where = ^{g^ — 2) and the numerals 8 and 3 in parentheses represent the statistical and 
systematic uncertainties in the last digits of the measured value. 1 ppm = 10~^. The world 
average value a^(exp) obtained from this and earlier measurements 0,0,0,0] is 

a^(exp) = 11 659 208 (6) x 10"^° (0.5 ppm). (2) 

This result provides the most stringent test of the Standard Model. 



Unfortunately, such a test must wait 



hadronic corrections to |6| 



'or further improvement in the uncertainty of the 



BSHHIiiHyHQ. The 



lowest-order hadronic 



vacuum-polarization effect has thus far been determined from two sources, (i) e^e~ annihi- 
lation cross section, and (ii) hadronic r decays. Several recent evaluations are listed in Table 
H] Their differences (except for the one obtained from the r decay data) are due to different 
interpretations and treatments of basically identical data. However, they all agree that the 
measurement of the e+e~ annihilation cross section, in particular in the region below p — uj 
resonances, must be improved substantially in order to reduce the experimental uncertainty 
significantly. Such efforts are underway at several laboratories. Particularly interesting and 
promising is new radiative-return measurements On the other hand, it is not clear at 
present whether the value from the r-decay data can be improved much further because of 
the difficulty in evaluating more precisely the effect of isospin breaking 0,0]- 

A new theoretical development is an attempt to calculate the hadronic vacuum- 
polarization effect on muon qf — 2 in lattice QCD 1181. 

nn 

The NLO hadronic contribution has been evaluated by two groups |8L I19I|: 
a^{had.NLO) = -10.1 (0.6) x 10~^°, 

a^ihad.NLO) = -9.8 (0.1)^.^ (0.0),,^ x lO'^". (3) 

The contribution from radiative corrections is identical in two papers. The small difference 
comes from the diagram in which two hadronic vacuum-polarizations are inserted in the 
second-order vertex diagram. 



TABLE I: Recent evaluations of lowest-order hadronic vacuum-polarization contribution to the 
muon g — 2. Some errors are separated according to their sources: measurement errors and radiative 
corrections. mentions a procedural error separately. 





process 


a^{had.LO) x 10^° 




Reference 


e" 


annihilation 


696.3 (6.2)e^,p(3.6)rad 




m 


e" 


annihilation 


694.8 (8.6) 




m 


e" 


annihilation 


692.4 {5.9)exp{2A)rad 






e" 


annihilation 


699.6 {8.5)e^p{l.9)rad 


(2.0)pj-oc 


M 


r 


decay 


711.0 (5.0)e..p(0.8),arf 


(2-8)sc/(2) 


m 



The contribution of hadronic light-by-light scattering to is more difficult to obtain a 
reliable value because it cannot utilize any experimental information and must rely solely 
on theory. After correction of a sign error, it seemed to have settled down to around 



on tneory. Alter cor 



a^(had.l - 1) ~ 80 (40) x 10"^^ (4) 
More recently, however, a considerably different value was reported \\^ : 

a^(had.l - 1) ~ 136 (25) x 10"", (5) 

which moves the prediction of the Standard Model closer to the experiment. This was 
obtained by imposing the short-distance QCD constraints on the '7r°7*7 amplitude, which 
was overlooked in previous analyses. Further confirmation of this result by a first principle 
calculation in lattice QCD would be highly desirable. 
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The weak interaction effect is known to two-loop order. The latest values are 

a^{weak) = 152 (1) x 10"" 

a^iweak) = 154 (1) (2) x 10"", (6) 
where (1) and (2) in the second line are the remaining theoretical uncertainty and Higgs 



mass uncertainty, respectively. Although the numerical difference between these values 
is insignificant for comparison with experiment, their approach to the fermionic triangle 
diagram seems to be different. We hope it is resolved before long. 

The QED contribution a^(QED), even though it is the predominant term of a^, has 
received little attention thus far because of its small error bars. The theoretical uncertainty 
comes predominantly from the term whose contribution to is about 3.3 ppm. The 
best value of a^(QED) reported previously (Eq. (11) of [2^) was 



a^(QED)o,rf = 116 584 705.7 (1.25)(1.15)(0.5) x 10"" 

= 116 584 705.7 (1.8) x 10"", (7) 

where 1.25 and 1.15 come from the uncertainties in the calculated and estimated ot' 
terms, respectively, and 0.5 is from the uncertainty in the fine structure constant a given in 
Eq. (17) of [22] obtained from the measurement and theory of Og. 

While updating ttf^^QED), however, we discovered that the previous evaluation of the 
term suffered from an error in a group of 18 Feynman diagrams This affects both Eq. 
(11) and Eq. (17) of Ref. [22] so that ((Tj) had to be revised. This discovery prompted us to 
reexamine all other terms contributing to a^j_{QED). 

The purpose of this paper is to report the result of this reexamination. We give a full 
account of 

(1) new evaluation of mass-dependent term of in an alternate formulation, 

(2) vastly improved numerical precision by an extensive numerical evaluation of 469 eighth- 
order Feynman diagrams, and 

(3) new evaluation of the term that depends on three masses me,m^,mT- (0.1094 (3) x 
10""), which replaces the old value (0.23 x 10~^^) quoted in \2M. 

n 

If one uses the latest value of a obtained from the atom interferometry measurement |2q| : 



a'^{a.i.) = 137.036 000 3 (10) [7.4 (8) 
the new estimate of the QED contribution becomes 

a^iQED) = 116 584 719.58 (0.02)(1.15)(0.85) x 10"", (9) 

where 0.02 replaces the previous uncertainty 1.25 of the term in ((Tj), an improvement of 
factor 60. The error 0.85 comes from the uncertainty in the fine structure constant Q;(a.i.) 
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given in (jSJ. Note that this error is larger than the corresponding error in ((Tj) because we 
used a;(a.i.) of (jHl) instead of the incorrect a{ae) used in ((Tj). The new value © is larger 
than (|7j) by 13.7 x 10^^^. Report on the improvement of Oe and a{ae) is being prepared j26|. 

As is seen from (0), the largest source of QED error is now the term, which was 
previously estimated to be 6.29(1.15) x 10^^^ 13 El- Although this is accurate enough for 
comparison with the current experimental data, a more precise value will become necessary 
in the future. It is being improved at present and will be reported shortly 28^ . 

Let us now present an outline of our approach to a^(QED) and a summary of results 
before going into details. The contribution of QED diagrams to can be written in the 
general form 



a^{QED) = Ai + A2{m^/me) + A2(m^/m^) + A3{mjme,mf,/mr), 



(10) 



where me, rrin, and rrir are the masses of the electron, muon, and tau, respectively. Through- 
out this article we shall use the values = 0.510 998 902(21)MeV/c2, = 105.658 
3568(52) MeV/c^, and = 1 777.05(29) MeV/c^, respectively |29|. 

The renormalizability of QED guarantees that Ai, A2, and A^ can be expanded in power 
series in a/vr with finite calculable coefficients: 



A?^ 



TT 



+ (- 



TT 



4.. 



1, 2, 3. 



:ii) 



A^^^ is known up to n = 4 from the study of the electron anomaly 2^|26|. AY\ ^^'^ 



l(2) /I (4) 



A^^ have been evaluated precisely by both numerical and analytic means. A^^^ is currently 



being improved by an extensive computer work [2a|. For the purpose of evaluating ci^( QED), 
however, we may use Ai obtained from the measured value of the electron anomaly 
subtracting small contributions due to muon, hadron, and weak interactions 



It is easy to see that A2 = A 



1(2) 



A 



(4) 



0: they have no corresponding Feynman 



diagram. A^\m^/mf,), A2^^(m^/me), and Af\m^/me,mn/mT-) have been evaluated by 
numerical integration, asymptotic expansion in mn/nie, power series expansion in me/m^, 
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and/or analytic integration. They are 



32 



33 



34 



Af\mf,/me,mij,/m, 



1.094 258 282 8 (98) 
7.8059 (25) X 10-^ 
22.868 379 36 (23), 
36.054 (21) X 10-^ 
52.763 (17) X 10-^ 



(12) 



where the errors are due to measurement uncertainty of and rrir only. The most striking 



feature of the term is the large size of Af\mn/mf.). It comes predominant 
diagrams involving a light-by- light scattering subdiagram, as was first discovered in 



from 
36l | and 

improved by numerical calculation [24]. Since Af\mn/me) is now known analytically js^, 
its uncertainty depends only on the uncertainty in the measurement of me/m^ and is totally 
negligible. 

The term A2 (m^/me) has been known by numerical integration only. A crude evaluation 
of contributing integrals made more than 10 years ago [2^, which was no more than an order 
of magnitude estimate, showed that A2 {m^/nie) contributes only about 3 ppm to a^. Thus 
it seemed that it was good enough for comparison with the experiment. Now that a program 
error was found in a part of evaluation of A2^''(m^/me) j2^ and since the measurement of 

is becoming more precise, however, it is important to re-examine these calculations and 
eliminate algebraic error, if any, completely and reduce the computational uncertainty as 
much as possible. 

Within the Feynman gauge two approaches had been developed for numerical integration 
of Feynman diagrams contributing to the anomalous magnetic moment [stI. An obvious 
and straightforward one is to evaluate each vertex individually and add them up. (This 
approach will be called Version B following j3].) Another one starts by combining several 
vertices into one with the help of the Ward-Takahashi identity 



g^A'^(p,g) = -S(p + |) + E(p-|) 



(13) 



where A'^{p,q) is the sum of vertices obtained by inserting the external magnetic field in 
fermion lines of a self-energy diagram S(p). p±q/2 is outgoing (incoming) muon momentum. 
Differentiating both sides of ()13p with respect to q,, one obtains 
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dqi 



Jg=0 

Obviously one may start from either the LHS or RHS of this equation to evaluate the 
anomalous magnetic moment. The approach based on the RHS and LHS will be called 
Version A and Version B, respectively. The former required some additional algebraic work 
but produced fewer integrals and ensured significant economy of computing time. 

Evaluation of the term was carried out in both Version A and Version i? jsTj. But, for 
the a"^ term, in particular for 126 diagrams containing a light-by-light scattering subdiagram, 
the Version B codes were so large that we chose initially to work only with Version A. For 
the reason discussed already we have now reevaluated them also in Version B We have 
now extended this effort to the remaining 108 diagrams and obtained their codes in Version 
B. Numerical evaluation shows that they are in good agreement with those of Version A. 
As a consequence all oc" diagrams contributing to a!^2^ {va^jme) have been confirmed by two 
or more independent formulations. We are confident that all codes are now free from any 
algebraic error. 

The remaining problem concerns the reliability of numerical integration. As a matter 
of fact, values of some integrals were called into question shortly after the old result was 
published 2^. It turned out that this was caused mainly by the relatively poor statistical 
sampling of the integrand resulting from shortage of computing power then available js^. 
The problem was made worse by the presence of severe non-statistical errors that originate 
from round-off errors inherent in all computer calculation. This will be called digit- deficiency 
errors. Various techniques had to be introduced to alleviate this problem {4^. See Appendix 
|B]for details. 

Now that the validity of codes is established we are justified to evaluate all integrals 
contributing to the term in either Version A or Version B, using vastly increased number 
of sampling points, made possible by the new generations of computers, and, at the same 
time, reducing digit- deficiency errors to a manageable level by various means. (See Appendix 

E) 

All integrals have been evaluated with successively increasing statistics over the period of 
more than 10 years. Some preliminary results were reported from time to time [4^. Only the 
latest and most accurate results are listed in Tables ITTl — IXHI Although earlier results are 
not shown explicitly, they have played crucial roles in checking the reliability of numerical 
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integration at every stages of calculation. 

The majority of integrals in the Version A calculation were found to be consistent with 
the results in Version B. But some of them were found to differ considerably because of poor 
statistical samplings and the d-d problem. Thus some Version A integrals have been reeval- 
uated to reduce the d-d problem. Evaluations of both versions are combined in quadrature, 
whenever appropriate, to improve the statistics. 

The latest value of yl2^''(m^/me) is 

APim^/me) = 132.682 3 (72), (15) 

which is larger by 5.2 than the old value 2^ 

A^^\m^/me) = 127.50 (41). (16) 

The difference between fjl5p and (fTBj) is partly accounted for by the correction of program 
error described in j^l but is mostly due to the fact that (fT^ suffered from poor statistics 
and the digit- deficiency problem. 

There is also a small contribution to from the three-mass term A^^\mn/me,m^/mj.) 
which arises from 102 diagrams containing two or three closed loops of v-p and/or /-/ type. 
Results of numerical evaluation are given in (plj) . (pT|) . and (jU^ . From these results we 
obtain 

APim^/me^m^/rrir) = 0.037 594 (83). (17) 



This is smaller than the value 0.079 (3) quoted in 2M- which corresponds to (jU^ obtained 
only from the diagrams containing /-/ loop, which were thought to be dominant. The new 
result (|17|) shows that this assumption was not fully justified. Another term of order is 
A^\mf^/mT-) which is calculable from 469 Feynman diagrams. However, its contribution to 

is of order {rrin/mrY ln{mr/mn)A^2\^) ~ 0.005 so that it may be safely ignored for now. 

Collecting all results of orders and js^l we find a^(QED) given in (jH)). In conclu- 
sion we have found that the improvement of the term does not significantly affect the 
comparison of theory and experiment of a^. The net effect of our calculation is to enhance 
the QED prediction ((7j) by 13.6 x 10^^^ and eliminate an important source of theoretical 
uncertainty. As far as QED is concerned, the term is now the most important source of 
uncertainty in a^. This is being improved The overall theoretical uncertainty of the 
Standard Model remains dominated by that of the hadronic vacuum-polarization effect. 
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II. CLASSIFICATION OF DIAGRAMS CONTRIBUTING TO Af\m^/me) 

There are altogether 469 Feynman diagrams contributing to Af{mjm^). Feynman 
integrals for these eighth-order vertex diagrams consist of twelve propagators integrated 
over 4 four- dimensional loop momenta. These diagrams have subdiagrams of vacuum- 
polarization {v-p) type and/or light-by- light scattering (/-/) type. The v-p subdiagrams 
found in /l2^''(m^/me) follows: 

112, which consists of one closed lepton loop of second-order. 

114, which consists of three proper closed lepton loops of fourth-order. 

114(2), which consists of three lepton loops of type 114 whose internal photon line has a 112 

insertion. 

rig, which consists of 15 proper closed lepton loops of sixth-order. 

The /-/ diagrams we need are: 
A4, which consists of six proper closed lepton loops of fourth-order, with four photon lines 
attached to them. 

A4 , which consists of 60 diagrams in which lepton lines and vertices of A4 are modified by 
second-order radiative corrections. 

We are now ready to classify the diagrams into four (gauge-invariant) groups: 

Group I. Second-order muon vertex diagrams containing lepton v-p loops 112, n4, 114(2) 
and/or He. This group consists of 49 diagrams. 

Group II. Fourth-order proper vertex diagrams containing lepton v-p loops 112 and/or 
114. This group consists of 90 diagrams. 

Group III. Sixth-order proper vertex diagrams containing a v-p loop 112. This group 
consists of 150 diagrams. 

Group IV. Muon vertex diagrams containing an /-/ subdiagram A4 with additional 2nd- 

(2) 

order radiative corrections, or one of A4 type. This group consists of 180 diagrams. 

All integrals of Groups I, II, and III have been evaluated by numerical means. Further- 



more, some of them have also been evaluated semi-analytically |43[. Group IV integrals 
have thus far been evaluated only by numerical integration, but in two independent ways, 
Version A and Version B, both in Feynman gauge. 

The starting point of Version A is the RHS of Eq. (fT^ . The algebraic structure of 
integrals in Version A is more complicated than that of Version B but their codes are 
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substantially smaller in general than the latter. For Group I, however, there is no advantage 
of using Version A. Thus this group is formulated in Version B only. 

All integrals are generated from a small number of templates, enabling us to make cross- 
checking of different diagrams, thereby reducing significantly possible programming errors. 
More information on Version A and Version B are given in Appendix lA 11 

Integrals thus obtained are divergent in general. Since computers are not capable of 
handling divergence directly, both ultraviolet (UV) and infrared (IR) divergences must be 
removed beforehand. We have introduced a two-step on-shell subtractive renormalization 
scheme, in which the first step removes both UV and IR divergences but does not give exact 
on-shell results. This is done to circumvent the inconvenient feature of the standard on-shell 
renormalization in which the renormalization terms do not remove and may even introduce 
extra IR-divergent terms. The second step yields the standard on-shell renormalization 
result when summed over all diagrams. 

The renormalization terms are generated in two ways: One by reduction of the original 
integral according to a well-defined power counting rule, and another from scratch, both 
analytically. This enables us to make extensive cross-checking between diagrams of various 
types and different orders. See Appendix lA 21 for more details. 

All integrals contributing to the term are evaluated numerically by the adaptive- 



44| . The major source of numerical 



iterative Monte-Carlo integration routine VEGAS 
uncertainty is the difficulty of accumulating a large number of good sampling points that 
do not suffer from the digit- deficiency problem caused by the round-off error. For this 
purpose quadruple precision is required in many cases. Unfortunately, this slows down the 
computation quite drastically. The accuracy of these integrals is checked by comparison 
with those obtained by other means whenever possible. The results of our calculation are 
summarized in the following sections. The reliability of these results, which depends critically 
on the reliability of the numerical integration routine VEGAS, is discussed in Appendix iBl 
Problems caused by non-statistical errors encountered in dealing with VEGAS and their 
solution are discussed there in detail. 



III. GROUP I DIAGRAMS 

Group I diagrams can be classified further into four gauge-invariant subgroups: 
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Subgroup 1(a). Diagrams obtained by inserting three 112 's (of electron/muon loop) in a 
second-order muon vertex. Seven Feynman diagrams belong to this subgroup. See Fig. ^a). 

Subgroup 1(b). Diagrams obtained by inserting a 112 and a 114 in a second-order muon 
vertex. Eighteen Feynman diagrams belong to this subgroup. See Fig. Wljo)- 

Subgroup 1(c). Diagrams containing 114(2). There are nine Feynman diagrams that belong 
to this subgroup. See Fig. |21 

Subgroup 1(d). Diagrams obtained by insertion of Hq in a second-order muon vertex. 
Fifteen Feynman diagrams belong to this subgroup. Eight are shown in Fig. El Diagrams 
a, c, (i, e, / and the time-reversed diagram of e have charge-conjugated counterparts. 

The evaluation of subgroups 1(a) and 1(b) is greatly facilitated by the analytic formulas 
available for the second- and fourth-order spectral representations of the renormalized photon 

n 

propagators [45]. The contribution to from the diagram obtained by sequential insertion 
of m k-ih order electron and n /-th order muon v-p loops into a second-order muon vertex 
is reduced to a simple formula 



ds 



Pk{s) 



4 l-y frue 
m 



dt 



Piit) 



4 l-y 



1 + 



, (18) 



where pk is the k-th order photon spectral function. Exact p2 and p4 can be found in Ref. 



38, 



47|. 



4^, 4^ . An exact spectral function for 114(2) and an approximate one for Hg are also available 




(b) 



FIG. 1: (a) Diagrams contributing to subgroup 1(a). (b) Diagrams contributing to subgroup 1(b). 
Solid horizontal lines represent the muon in external magnetic field. Numerals "2" , "4" within solid 
circles refer to the proper renormalized v-p diagrams 112 and 114 ! respectively. Letters li, hjh refer 
to electron or muon. Seven and 18 Feynman diagrams contribute to 1(a) and 1(b), respectively. 

The contribution of diagrams of Fig. [T]can be obtained by choosing {k = 2,m = 3,n = 

11 



0), (fc = 2, m = 2, / = 2, n = 1), (/c = 2,m = 1,/ = 2,n = 2). The latest numerical values 



obtained by evaluating these integrals using VEGAS are listed in Table |nj where the 
number of sampling points per iteration and the number of iterations are also listed. 

Note that these diagrams need no additional renormalization. Thus the renormalized 
amplitudes a^vi^li given by 



Note also that, to be consistent with notations used later, M^^^2^^ , etc., should have been 
written as M2^^2.^'^\ etc. The first superscript /i is often (but not always) suppressed for 
simplicity when there is no danger of confusion. 

Adding up the first three rows of Table m we obtain the total contribution of diagrams 
of subgroup 1(a) 

ai'l, = 7.745 140 (30) . (20) 

This is abou. 40 «,„e. .o.e p.eCe .he ea... .e.u. Q. Fu.he™o.e iti, i„ excellent 

agreement with the results obtained by an asymptotic expansion in m^/mg 



43|: 



afl^{asymp) = 7.745 136 8 (8), (21) 

where the uncertainty comes only from the measurement of muon mass. 

The contributions of Fig. [Tfb) for {li,^) = (e, e), (e, /x), and (/i, e) can be written down 
in a similar fashion. The most recent results of numerical integration by VEGAS are listed 
in the last three rows of Table IHl These diagrams need no additional renormalization, too. 
The sum of these results is the contribution of the subgroup 1(b) 

af^^ = 7.581 262 (50). (22) 

This again is in excellent agreement with the asymptotic expansion result 

4^[)(as?/mp) = 7.581 275 5 (2), (23) 

where the uncertainty comes only from the muon mass. 

In evaluating the contribution to from the 9 Feynman diagrams of subgroup 1(c) shown 
in FIG. 121 our initial approach was to make use of the parametric integral representation 

(2) 

of the v-p term 114 . Following the two-step renormalization procedure, these contributions 
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TABLE II: Contributions of diagrams of Figs. E^a) and E^b). np is the number of Feynman 
diagrams represented by the integral. These evaluations were carried out on a workstations in 
2001. 



Integral np Value (Error) Sampling per No. of 

including np iteration iterations 

M^y^':l 1 7.223 077 (29) 1 x 10^ 100 

^2!p2:3 3 0.494 075 ( 6) 1 x 10^ 100 

^2rP2':3 3 0.027 988 ( 1) 1 x 10® 60 

^iylpA 6 7.127 996 (49) 1 x 10^ 100 

^2!p2.pa 6 0.119 601 ( 3) 1 X 10® 60 

^iylpA 6 0.333 665 ( 4) 1 x 10® 60 





FIG. 2: Diagrams contributing to subgroup 1(c). (hjh) = (e, e), (e, /i), or {fi,e). See FIG. 1 for 
notation. 



can be written in the form 14811. 



^(8) _ Jluh) 



^2,P4(P2)' (2'^) 

('1,^2) 



V(C) 

where each term of 

"2,P4(P2) ~ ^^'^2,P4a(P2) ^ '^^^'^2,P46(P2) ^^-'-'2,P2 ^'^2,P2 ; ) 

are finite integrals obtained by the K5 renormalization procedure described in Ref. and 
Appendix A. The suffix P2 stands for the second-order v-p diagram 112, -P4 for the fourth- 
order v-p diagram 114, while P4(P2) represents the diagram 114(2) • -P4 receives contributions 
from P4a (vertex correction) and (lepton self-energy insertion), P4 = P4a + 2P4b. 
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TABLE III: Contributions of diagrams of Fig. [21 np is the number of Feynman diagrams repre- 
sented by the integral. Numerical work was carried out on a workstations during 2001. 



Integral np Value (Error) Sampling per No. of 



2,P4a(P2) 

2,P4a{P2) 
f(e,e) 
^2,P46(P2) 

^2,P46(P2) 

^2,P46(P2) 



AMi't;':L^.^ 2 



including rii? 


iteration 


iterations 


0.597 477 1 (111) 


1 X 10^ 


100 


0.121 902 1 (58) 


1 X 10"^ 


100 


0.021 017 1 (13) 


1 X lO'^ 


100 


0.982 017 4 (109) 


1 X 10^ 


100 


0.099 244 1 (84) 


1 X lO'^ 


100 


0.000 586 ( 4) 


1 X lO'^ 


100 



The results of numerical evaluation of ()25p . obtained by VEGAS, are listed in Table 
nil Numerical values of lower-order Feynman integrals, in terms of which the residual 
renormalization terms are expressed, are given in Table IIVI From these tables and ()24|1 we 
obtain 

4>1(P2) = 1-440 744 (16), (26) 
4!pi(P2) = 0-161 982 (11), (27) 



«K(P2) = 0-021 583 (2). (28) 

The new results fl26|) and ()28|) confirm the old results but with a much higher precision. For 
(jTTj) the agreement between the old and new values is rather poor. 

About a decade ago the leading log term of ci^2'pi{P2) obtained by the renormalization 

Hup method {3] seemed to disagree with the numerical evaluation. However, it was found 
that this was caused by an improper use of the asymptotic photon propagator obtained 
for massless QED in It is important to note that the asymptotic photon propagator 
for massless QED is not the same as one for massive QED as was proven explicitly in Isoj ]. 
Use of the correct photon propagator in the renormalization group method leads to results 
which agree very well with the numerical integration result (52l. l53l| . This episode provides 
an explicit example of danger of confusing the asymptotic behavior with the mass-less limit, 
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TABLE IV: Auxiliary integrals for Group I. Some integrals are known exactly. Some are obtained 
by expansion in me/m^ to sufficiently high orders. Their uncertainties come from that of mg/m^ 
only. The remaining integrals are obtained numerically by VEGAS. Total sampling points are of 
order 10^^. 

Integral Value (Error) Integral Value (Error) 

M^^^j^l 1.094 258 282 7 (98) M^^^^^ 0.015 687 421 •• • 

^2^P2* -0.161 084 05 ••• 

AS2 0.75 ^^fra 0-063 399 266 •• • 

AS^I'^ 1.885 732 6 (158) A-B2!p2 9-405 5 x 10"^ 

AL4 0.465 024 (12) -0.437 094 (21) 

A(5m4 1.906 340 (22) 



which results in different non-leading terms. 

We obtained an independent check of ()26|) using an exact spectral function for 114(2) of 



Fig. 121 which was derived 5J] from the QCD spectral function obtained in 



3- 



Numerical 



integration using this spectral function gives 

^ip\{P2) = 1-440 622 (173), (29) 

for 100 million sampling points iterated 100 times in quadruple precision. This is in agree- 
ment with (j26|l to the fifth decimal point although their approaches are completely different. 
Undoubtedly both and must be correct. 

The best value of af)s is obtained by adding up (pUj) . (j27j), and (PHj) : 



/(c) 



a 



1.624 308 (19). (30) 



The contribution to from 15 diagrams of subgroup 1(d) (see Fig. Oj) can be written as 

a2,P6i = AM2,p6i + residual renormalization terms , {i = a, . . . , h). (31) 

Divergence- free integrals AM2,p6j are defined by (4.13) of Ref. j^l- Their numerical values 
(summed over the diagrams related by time-reversal and charge-conjugation symmetries) 
are evaluated numerically by VEGAS and listed in the third column of Table 

Summing up the contributions of diagrams a to /i of Fig. El we obtain the following 
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P6a P6b P6c P6d 




P6e P6f P6g P6h 



FIG. 3: Eighth-order vertices of subgroup 1(d) obtained by insertion of sixth-order (single electron 
loop) v-p diagram Xlg in a second-order muon vertex. 

TABLE V: Contributions of diagrams of Fig. El rii? is the number of Feynman diagrams repre- 
sented by the integral. M2^pQe was evaluated in 1998 on SP2 at Cornell Theory Center. Others 
were evaluated in 1998 on Fujitsu VX at Nara Women's University, Japan. 

Integral np Value (Error) Sampling per No. of 







including np 


iteration 


iterations 


AM2,P6a 


2 


5.676 002 (168) 


4 X 10^ 


60 




1 


3.058 301 (152) 


2 X 10^ 


60 


AM2,P6c 


2 


1.483 501 (104) 


2 X 10^ 


60 


AM2,P6d 


2 


-3.127 282 (122) 


2 X 10^ 


60 


AM2,P6e 


4 


-0.073 885 (234) 


6 X 10^ 


60 


AM2,P6/ 


2 


-4.064 113 (151) 


2 X 10^ 


60 


AM2,P6g 


1 


-0.247 237 (100) 


2 X 10^ 


60 


AM2,P6h 


1 


2.838 657 ( 74) 


2 X 10^ 


60 



expression: 



4w = EL^.AM2,p6, - 4AB,AM^^i 

- 2A6m4M^%%, (32) 
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where 



A'B, 



AM 



2,P4 
AL4 

A(5m4 



AM 



+ 2AM 



3 

4' 



'2,P4a ^ ^^^''-'2,P4b' 

AL4^ + 2AL4C + AL4i + 2AL4S, 

= A_B4a + AB^h, 

= Adm^a + ASrriib- 



(33) 



The quantities hsted in are defined in Ref. 40]. Their numerical values are listed in 



Table HVl The 1998 results of numerical integration of AM2,p6i are listed in Table 
the numerical values in Tables |TV| and |V| we obtain the value reported previously 



From 



h{d) 



0.230 596 (416). 



(34) 



This deviates strongly from the old result —0.7945(202) 2J|. The problem with 



was first pointed out in jssi in which a^^^^ was evaluated without the 0{m^/m^) term by 

a renormalization group method. Soon afterwards a Fade approximant of the sixth-order 

I I 

photon spectral function was used to evaluate the full correction |55l] : 



{Fade) 



0.230 362 (5). 



(35) 



Our new result (jH^ is in good agreement with The primary cause of the old discrepancy 
was traced to very poor statistics of the original evaluation 15611 . Increase of statistics by 
two orders of magnitude improved the result to —0.2415(19) [39]. However, the discrepancy 
with ()35p was still non-negligible. Finally, the problem was traced to round-off errors caused 
by insufficient number of effective digits in real*8 arithmetic in carrying out renormalization 
by numerical means This was resolved by going over to the real*16 arithmetic. (See 
Appendix IHl for further discussion on this point.) 

Note that the uncertainty in (|35|) may be an underestimate since it does not include 
the uncertainty of the Fade approximation itself. However, it seems to be small compared 
with the quoted uncertainty 40|. In principle it is possible to prove or disprove it by more 
numerical work. However, it would require 6,000 times more computing time in order to 
match the precision of achieved by the Fade method. This is not only impractical but 
also pointless since there is no need to improve the current precision further. 
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Collecting the results (EH), (P^. (IHUj) and we find the best value of the contribution 
to the muon anomaly from the 49 diagrams of group I: 

a^^j^ = 16.720 359 (20) . (36) 

IV. GROUP II DIAGRAMS 

Diagrams of this group are generated by inserting 112 and 114 in the photon lines of fourth- 
order muon vertex diagrams. Use of analytic expressions for the second- and fourth-order 
spectral functions for the photon propagators and time-reversal symmetry cuts down the 
number of independent integrals in Version A from 90 to 11. 

The contribution to arising from the set of vertex diagrams represented by the "self- 
energy" diagrams of Fig. |3] can be written in the form 

a^^P^ = AM'4 + residual renormalization terms, (37) 

where AM4 are finite integrals obtained in the intermediate step of two-step renormaliza- 
1 I 

tion 57]. Their numerical values, obtained by VEGAS are listed in Table IVIl The values of 
auxiliary integrals needed to calculate the total contribution of group II diagrams are given 
in Tables HVl and IVTTl 




FIG. 4: Eighth-order diagrams obtained from the fourth-order vertex diagrams by inserting 
vacuum-polarization loops 112 and 114, which consist of either electron or muon loop. 

Summing the contributions of diagrams of the first, second, and third rows of Fig. HI one 
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obtains 



a4,P4 = 2AMij;U AM^^J)^,^ + AM^'^^ 



4f),Pl':4 ~'~ ^^^^4b,P0:4 

- AB^M^p'l - AB^ylM2, (38) 

«4.^2,P2 = AMi:;$2,P2 + ^<:pi'.2,P0:2 

_ A r(^''^) A/f ('^'^^ 
'-^^2,P2 ^^-'2,P2 

-h ZZAJW4„ p2,P2 + ^^'^46,P1':2,P0:2 + ^^^-'4b,Pl' :2,P0:2 

- AE^:^'^^^ M!i:,i - ABi^i Mif/,\ (39) 

a4,P2:2 = 2AMt%.^, + AM[1';'X':2:2 + ^<P0:2:2 

- AE^M^J,!^ - A^fpL ^2 

+ 4AMi:;^)2:2 + 2^M^Zpi':2:2 + 2 A ^1,^;^;)^ 2:2 

- 2Ai?2M(;p';)2 - 2^Biyi, M2, (40) 

respectively, where M2 is equal to AM^fj^} -2AB2M!fj^^ . The factor 2 in front of AM^a,... 
accounts for equivalent diagrams obtained by time-reversal and another factor 2 in front of 
AMia,... and AM4b^... accounts for interchange of electron and muon vacuum-polarization 
loops. In contrast, the auxiliary integrals listed in Tables HVl and Ivnl do not include multi- 
plicity. Following the convention adopted below Eq. ()19p . the first superscript /i indicating 
the external muon line is supressed for simplicity. For instance, AM^'^pl^ p2 is written as 

^^'^4a,P2,P2- 

Substituting the data from Tables HYI and EH into dHHl), dHHI), and (001) we obtain 

a4_P4 = - 2.778 565 (253), 
a4,P2,P2 = -4.553 017 ( 68), 
a4,P2:2 = -9.342 599 (438). (41) 

Two of these terms were also evaluated in |i4^] by an asymptotic expansion in m^/mg: 

a4,P4{asymp) = - 2.778 852 33 (5), 
a4,P2:2{asymp) = -9.342 722 1 (5). (42) 

They are in excellent agreement with the numerical integration results. 
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TABLE VI: Contributions of diagrams of Fig. 0] np' is the number of Feynman diagrams repre- 
sented by the integral, t.r. refers to time-reversed amplitude. Numerical evaluation was carried 
out on a workstations in 2001. 



Integral 


np 


Value (Error) 


Sampling per 


No. of 






including np 


iteration 


iterations 


^M^X + t.r. 


18 


2.047 838 (221) 


1 X 10^ 


100 














18 


-2.486 595 (119) 


1 X 10^ 


120 




3 


2.289 959 (144) 


1 X 10^ 


100 




6 


0.054 120 ( 34) 


1 X 10^ 


100 


A )i#(e,e) 
^^'^4b,Pl':2,P0:2 


3 


-4.249 598 ( 76) 


1 X 10^ 


100 


A 

^^"46,P1':2,P0:2 










T^^^"4fe,Pl':2,P0:2 


6 


-0.485 108 ( 14) 


1 X 10^ 


100 


AMi:f;2:2 + t-r. 


6 


5.148 441 (377) 


1 X 10^ 


100 




12 


0.260 977 (103) 


1 X 10^ 


100 


^-'*^46,P0:2:2 










T^^^^^46,P1':2:2 


6 


-8.633 608 (190) 


1 X 10^ 


100 


^-'^-'46,P0:2:2 










+ AMi,^;^),2:2 


12 


-1.102 819 ( 42) 


1 X 10^ 


100 



TABLE VII: Auxiliary integrals for Group II. Some integrals are known exactly. Some are obtained 
by expansion in rrielvn^ to necessary orders. Their uncertainties come from that of rrie/m^ only. 
Remaining integrals are obtained by VEGAS integration, with total sampling points of order 10^^. 

Integral Value (Error) Integral Value (Error) 

Ma 0.5 M^^^l 1.493 671 581 (8) 

^2!p2:2 2.718 655 7 (1) M^^^j^l^ 0.050 259 648 (1) 

A52 0.75 ^B^yi 0.063 399 266 .. . 

AB^yl^ 5.330 381 (61) ^B^'^f^.^ 0.236 018 (9) 

AS^''/^^ 2.439 109 (53) 
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Combining the results (jl^ and the value of a4,P2,p2 from (PT|) we find the best value for 
the contribution of 90 diagrams of group II to be 

= - 16.674 591 ( 68). (43) 



V. GROUP III DIAGRAMS 



Diagrams belonging to this group are generated by inserting a second-order vacuum- 
polarization loop 112 in the photon lines of sixth-order muon vertex diagrams of the three- 
photon-exchange type. Time-reversal invariance and use of the function p2 (see (jl8j) ) for the 
photon spectral function reduce the number of independent integrals in Version A from 150 
to 8. Some of these integrals are represented by the "self-energy" diagrams of Fig. El 




FIG. 5: Typical eighth-order diagrams obtained by insertion of a vacuum-polarization loop 112 in 
muon diagrams of the three-photon-exchange type. Altogether there are 150 diagrams of this type. 

Let Mga^p be the magnetic moment projection in Version A of the set of 150 diagrams 
generated from a self-energy diagram a (=A through H) of Fig. El by insertion of 112 and an 
external vertex. The renormalized contribution due to the group III diagrams can then be 
written as 

H 

^ni = '^Vaafia,P2, (44) 

a=A 

where 

0'6a,P2 = A^6a,P2 + rcsidual renormalization terms. (45) 



where all divergences have been projected out by K5 and Iji operations. (See Ref. [57[.) 

The latest numerical values of Group III integrals are summarized in Table IVlTTl Numer- 
ical values of auxiliary integrals needed in the renormalization scheme are listed in Tables 
IIV| I VI II and |TXl For comparison, the results of old calculation [2^ carried out in dou- 
ble precision are listed in the last column of Table IVIIII This is to examine the effect of 
digit-deficiency error. In this case the effect is relatively mild because the introduction of a 
vacuum-polarization loop tends to make the integrand less sensitive to the singularity. 
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TABLE VIII: Contributions of diagrams of Fig. \^ np is the number of Feynman diagrams 
represented by the integral. This calculation was carried out in quadruple precision in 2001 
2003 on a workstations to examine the influence of digit-deficiency error in the calculation of 
carried out in double precision. 



24| 



Integral 


np 


Value (Error) 


Sampling per 


INo. OI 


Data from 






including np 


iteration 


iterations 


Ref. [24j 


AMQa,P2 


15 


-12.934 780 (1081) 


4 X 10^ 


100 


-12.940 1 (130) 




15 


18.797 294 (1309) 


4 X 10^ 


140 


18.797 (171) 


AM6c,P2 


15 


3.997 996 (1773) 


4 X 10^ 


100 


4.000 7 (178) 


AMed,P2 


30 


10.492 627 (1507) 


8 X 10^ 


111 


10.494 (225) 


AM6e,P2 


15 


10.990 435 ( 981) 


4 X 10^ 


119 


11.000 1 (121) 


AM6/,P2 


15 


5.652 451 (1503) 


4 X 10^ 


100 


5.651 8 (166) 


AM6g,P2 


30 


19.747 805 (1558) 


4 X 10^ 


100 


19.742 4 (172) 


AMeh,P2 


15 


-18.363 491 (1433) 


4 X 10^ 


100 


-18.361 5 (141) 



When summed over all the diagrams of group III, the UV- and IR-divergent pieces cancel 
out and the total contribution to can be written as a sum of finite pieces: 



(8) 



+ {M^2^%[I] - M^rj2)A5m^ + (M24/] - M2^)A5m[^^l 
- M^^^j^l[ABi + 2ALi - 2(A52)2] 
- M2{AB[^f^+2AL^^f^ - AAB^AB^^f^). (46) 

Plugging the values listed in Tables llV | IVIII and IIXI in (j46|l . we obtain 

af}j = 10.793 43 (414). (47) 

The error in ()47|) can be reduced easily if necessary. The ration a!f]j/cS^\ where a(6) is the 
value of sixth-order muon moment without closed lepton loop, is about 11, which is not 
very far from the very crude expectation "iK ~ 9, where K is from (jHH|) . although such a 
comparison is more appropriate for individual terms on Table IVTTTI than their sum. See Sec. 
VIII for further discussion of enhancement factor K. 
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TABLE IX: Auxiliary integrals for Group III. Some integrals are known exactly. Some are obtained 
by expansion in mg/rra^ to necessary orders. Their uncertainties come from that of rrie/m^ only. 
Remaining integrals are obtained by VEGAS integration, with total sampling points of order 10^^. 



Integral Value (Error) Integral Value (Error) 

M2* 1.0 M2.[/] -1.0 

M2^^p'2 2.349 621 ( 35) M^^'^)^[I] -2.183 159 (95) 

AM4 0.030 833 612 ... AAfj^/g^ -0.628 831 80 (2) 

l^Lfp2 3.118 868 (201) ^sf 'pl -3.427 615 (237) 

A(5m4 1.906 340 (22) AJm^^],^^ 11.151 387 (303) 



VI. GROUP IV DIAGRAMS 



Diagrams of this group can be divided into four subgroups: IV(a), IV(b), IV(c), and 
IV(d). Each subgroup consists of two equivalent sets of diagrams related by charge conjuga- 
tion (reversal of the direction of momentum flow in the loop of the light-by- light scattering 
subdiagram). Diagrams of subgroups IV(a), IV(b), and IV(c) are obtained by modifying 
the sixth-order diagram which contains the light-by- light scattering subdiagram A4, one of 
whose external photon line represents the magnetic field. The magnetic moment contribu- 
tion Mqll of this sixth-order diagram is known analytically 34], whose numerical value is 

MgLL = 20.947 924 34 (21) (48) 

when A4 is an electron loop, and the uncertainty is due to that of the muon mass only. 

Subgroup IV(a). Diagrams obtained by inserting a second-order vacuum-polarization 
loop 112 in Mqll. They are all appropriate modifications of the integral Mqll,P2 defined 
by (2.4) of Ref. js^. Denote these integrals as Afg^^pg where {hjh) = (e, e), (e,/i) or 
(/i,e). This subgroup is comprised of 54 diagrams. They are generically represented by the 
self-energy-like diagrams shown in Fig. IHl 

Subgroup IV(b). Diagrams containing sixth-order light-by- light scattering subdiagrams 
Ag. Altogether, there are 60 diagrams of this type. Charge-conjugation and time- reversal 
symmetries and summation over external vertex insertions reduce the number of independent 
integrals to 4 in Version A. These integrals are generically represented by the self-energy-like 
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(a) 



(b) 



(c) 



FIG. 6: Muon self-energy- like diagrams representing the external- vertex-summed integrals of 
subgroup IV(a). {li,^) = (e, e), (e, /x), or e), where li, I2 refer to the light-by-light scattering 
loop A4 and vacuum-polarization loop 112, respectively. 





LLL 



FIG. 7: Muon self-energy- like diagrams representing (external- vertex-summed) integrals of sub- 
group IV(b), IV(c), and IV(d). 

diagrams LLA, LLB, LLC and LLD of Fig. [3 

Subgroup IV(c). Diagrams obtained by including second-order radiative corrections on 
the muon line of Mql^. There are 48 diagrams that belong to this subgroup. Summation over 
external vertex insertions and use of the interrelations available due to charge-conjugation 
and time-reversal symmetries leave five independent integrals in Version A. They are gener- 
ically represented by the self-energy-like diagrams LLE, LLF, LLG, LLH and LLI of Fig. 

m 
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(49) 



Subgroup IV(d). Diagrams generated by inserting A4 internally in fourth-order vertex 
diagrams. Diagrams of this type appear for the first time in the eighth order. Charge- 
conjugation invariance and summation over the external vertex insertion with the help of 
the Ward-Takahashi identity lead us to three independent integrals in Version A. They are 
represented by the diagrams LLJ, LLK and LLL of Fig. [7| No further discussion of this 
subgroup will be given in this paper since it was treated in a separate paper [2^. 

In subgroups IV(a), IV(b), and IV(c) UV-divergences arising from the light-by-light scat- 
tering subdiagram A4, or more explicitly W^l^^iq, ki, kj, k), can be taken care of by making 
use of the identity: 

'_d_ 

which follows from the Ward-Takahashi identity. Namely, no explicit UV renormalization is 
needed if one uses the RHS of fl49|l instead of LHS and the fact that of (fT^ vanishes 
by Furry's theorem. On the other hand, is nonzero for subgroup IV(d) and the UV- 

divergence associated with the light-by- light scattering subdiagram A4 must be regularized, 
e. g., by dimensional regularization. For these diagrams it is necessary to carry out explicit 
renormalization of A4 as well as that of the two sixth-order vertex subdiagrams containing A4. 
See j^l for a detailed discussion of renormalization based on a combination of dimensional 
regularization and Pauli-Villars regularization. 

As was announced in Sec. all diagrams of Group IV have now been evaluated in 
both Version A and Version B. In the following let us consider Version A and Version B 
separately since renormalization is handled slightly differently in two cases. 

A. Version A 

The calculation of group IV(a) contribution is particularly simple. This is because Mq^l 
has been fully tested by comparison with the analytic result 2J|, and insertion of the vacuum 
polarization term is straightforward. One therefore finds that integrals Mg^]| are all finite, 
which means 
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TABLE X: Contributions of diagrams of Fig. [HI np' is the number of Feynman diagrams represented 
by the integral. The main term AMg2'2''p2 evaluated on vl at Cornell Theory Center in 2001. 
The rest was evaluated in 2001 on Condor cluster at University of Wisconsin. 



Integral 


np 


Value (Error) 


Sampling per 


No. of 






including np 


iteration 


iterations 


A /r/-(e.e) 


18 


116.759 183 ( 292) 


6 X 10^° 


180 


A /i^(e>/^) 


18 


2.697 443 ( 142) 


1 X 10^ 


110 


A /i#(M.e) 


18 


4.328 885 ( 293) 


1 X 10^ 


100 



where /i, I2 refer to the light-by-light scattering loop A4 and vacuum-polarization loop 112, 
respectively. Thus the contribution of subgroup IV(a) can be written as 

= E ^<tp2' (51) 

where the individual terms are given in Table |Xl 

Let us denote magnetic projections of subgroups IV(b) and IV(c) as Msna where a = 
A, ..,1. Relating the IR- and UV-divergent Mgna to the finite, numerically calculable piece 
AMsLLa defined by the procedure of two-step renormalization of Ref. j^, one can write 
the contributions of the diagrams of subgroups IV(b) and IV(c) as 

D 

"iVib) ~ 



afL. = J2r]^AMsLLa - SAB^M^ll, (52) 



a=A 

and 

4v(c)= I^^aAMsz^ia - 2A52M6ii. (53) 

a=E 

JNTumerical integration of all terms contributing to a^y has been carried out using VEGAS 



4J]. The latest results for Groups IV(b) and IV(c) are listed in Table IXll The result for 
Group IV(d) had been handled separately Qj. In general, the major difficulty in dealing 
with the diagrams of Groups IV(b) and IV(c) arises from the enormous size of integrands 
(up to 5000 terms and 240 kilobytes of FORTRAN source code per integral) and the large 
number of integration variables (up to 10). 

Diagrams of Groups IV(b) and IV(c) have singular surfaces just outside of the integration 
domain (unit cube) at a distance of ~ (mg/m^)^ ~ 1/40, 000. This makes the evaluation of 
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TABLE XI: Contributions of diagrams of Fig. Qexcluding LLJ, LLK, and LLL which were evaluated 
separately in Ref. 



23l |. np is the number of Feynman diagrams represented by the integral. Some 
integrals are split into two parts: d-part is evaluated in real*8 and q-part is evaluated in real*16. 



42| . The superscript * indicates 



a-part refers to the adjustable precision method developed by 
that indicated contributions were obtained by extrapolation from calculations in which the edges 
of integration domain were chopped off by l.d-10. See Appendix^for details. Numerical work was 
carried out on SP3, velocity cluster, SP2, Condor cluster, and a workstations over several years. 
The table lists only the latest of results obtained by various means. 



Integral 


np 


Value (Error) 


Sampling per 


No. of 






including np 


iteration 


iterations 


AMsLLA 


10 


52.063 459 (1497) 


4 X 10® 


3300 


AMsLLB 


20 


- 75.014 508 (1838) 






d-part 




- 53.000 600 ( 981) 


1 X 10^° 


430 


a-part* 




- 22.013 908 (1554) 


4 X lO'^ 


460 


^MsLLC 


20 


107.488 810 (2811) 


4 X 10® 


5900 


^MsLLD 


10 


- 37.824 352 (1137) 


1 X 10^° 


200 


AMsLLE 


6 


- 21.607 656 (1053) 






d-part 




- 20.920 745 ( 446) 


1 X 10^° 


304 


a-part 




- 0.686 911 ( 954) 


2 X lO'^ 


280 


AMsLLF 


12 


- 75.765 816 (2341) 


1 X 10^° 


1000 


AMsLLG 


12 


- 35.077 389 (1410) 


1 X 10^° 


470 


AMsLLH 


6 


54.025 704 (2411) 






d-part 




51.820 951 (889) 


2 X 10^° 


470 


q-part* 




2.204 753 (2241) 


4 X 10"^ 


391 


AMsLLI 


12 


112.756 785 (2683) 


1 X 10^° 


450 



their contributions to A2 ijf^n/t^e) much more sensitive to the d-d problem compared with 
the evaluation of the same set of diagrams contributing to the mass-independent A\ whose 
singularity is far outside (~ 1) of the domain of integration. 

Because of the d-d problem intensified by the proximity of the singularity, all strategies 
discussed in Appendix |Bl had to be tried to evaluate these integrals. 

In most cases the first step is to make the integrand smoother by stretching (see Appendix 
IB|) . which is repeated several times until the integrand behaves more gently. 
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Although chopping (see Appendix |B)) was handy to obtain a rough estimate quickly, we 
had to abandon it in the end because extrapolation to 5 = turned out to be too unreliable 
in order to reach the desired precision. 

Most integrals were then evaluated by splitting them into two parts, one evaluated in 
real*8 and the other in real*16. In some cases, however, even the part evaluated in real*16 
suffered from severe d-d problem, preventing us from collecting large enough samplings for 
high statistics. Analyzing this problem closely, we found that it is possible to evaluate 
these integrals by the following procedure: First try several iterations with a positive rescale 
parameter f3 (typically f3 = 0.5) until VEGAS begins to show strong sign of blowing up 
due to the d-d problem. Then freeze (3 to (see Appendix^. This may solve the problem 
in most cases. If not, try several iterations and see how rapidly the calculation runs into 
the d-d problem. It turns out that it takes place very early if we chose too many sampling 
points Ns per iteration. This is because choosing large Ns increases the chance of hitting 
random numbers too close to the singularity within one iteration. As a consequence the 
d-d problem is likely to dominate each iteration and makes it very difficult to collect large 
enough number of good samplings. We found that a better strategy is to reduce the size 
of Ns to a moderate value and, instead, increase the number of iterations Nj substantially. 
This is acceptable since, for (3 = which means that the distribution function p is no longer 
changed from iteration to iteration, the final error generated by VEGAS depends only on 
the product NsNj. 

This strategy has been applied in particular to the diagrams LLA and LLC, as is seen 
from Table IXll Entries in Table IXll are only the best of results obtained by various methods 
discussed above. They are consistent with each other despite their diverse approaches. 

One obtains from Tables 1X1 and IXll the contributions of subgroup IV(a) and the Version 
A contributions of subgroups IV(b) and IV(c): 



a 



(8) 

'IV{a) 



123.785 51 ( 44), 



a 



(8) 

'IVib) 



0.419 42 (385), 



a 



(8) 

'IV{c) 



2.909 41 (459). 



(54) 
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B. Version B 



In Version B the magnetic moment projection is evaluated for each vertex diagram on the 
LHS of (jl4j) . It is convenient to denote these diagrams in terms of self-energy-like diagrams 
of Fig. [7[ by attaching suffix i to indicate the lepton line in which an external magnetic field 
vertex is inserted. For instance, we obtain vertex diagrams LLAl, LLA2, LLA5 from the 
diagram LLA. 

We will not discuss subgroup IV(a) here since its Version A has already been fully tested. 
For subgroup IV(b) we find 

D 5 

= 5^5^r/,AM8ii,, - A^B^M^LL. (55) 

a=A 1=1 

instead of ()52|). Note that the last term of ()55p is different from that of (jH2I)- This is not an 
error. It arises from difference in the definition of AM terms. 
Similarly, for subgroup IV(c) we obtain 

/ 3 

4?(e)= ^^r/aAMsLLai - 2A52M6LL. (56) 

a=E i=l 

The results of numerical evaluation are listed in Table IXIII Precision of these calculations 
is still modest but high enough to show the consistency with the calculation of Version A. 
See the last column of Table IXIII for comparison of two Versions. Numerical work has been 
carried out with the same care as that described for Version A. The numerical calculation 
of AMsLLB was particularly difficult. 

One obtains from Table IXTTl the values of afy^^^^ and ct/y^,,) 

= - 0.372 (168), 
af^(^) = 2.876 3 (173), (57) 

which are not inconsistent with those given in (jHH), although much improvement is needed 
to become competitive with the Version A results. 

C. Total contribution of Group IV 

The contribution of subgroup IV(a) is listed only in ()54|) since it was not evaluated in 
Version B. The statistical combination of two versions of subgroups IV(b) and IV(c) is 
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TABLE XII: Contribution of Group IV(b) and Group IV(c) diagrams of Fig. [TJevaluated in Version 
B. Double precision is used for all calculations which are carried out on Fujitsu VPP at RIKEN. 
Finite renormalization terms AB2MQLLi,i = 1,2,3, are needed for LLA and LLC, respectively, in 
order to compare them with the calculations in Version A. Mqii(^2+3) = -^^6LL2 + AfeLLS is obtained 
subtracting Mqlh from the known value of Mqll{= Mqlli + ^6LL2 + -^eLLs) given in 



Integral 




Value (Error) 


Sampling per 


No. of 


Difference 






including np 


iteration 


iterations 


Ver.A - Ver.B 




10 


52.080 79 ( 731) 






-0.016 18 ( 785) 


YJl^i r]AAM8LLAi 




60.467 98 ( 731) 


1 X 10^ 


122 




—AB2MQLL2 




-8.387 20 ( 15) 


1 X 10^° 


280 




^MsLLB 


20 


-74.999 66 (1060) 


1 X 10^ 


544 


-0.014 83 (1076) 


AMsLLC 


20 


107.503 69 ( 877) 






-0.012 44 ( 955) 


ELi VC^MsLLCi 




114.827 43 ( 876) 


1 X 10^ 


357 




-Ai?2M62.L(l+3) 




-7.323 75 ( 15) 








AMsLLD 


10 


-37.823 98 ( 580) 


1 X 10^ 


120 


-0.000 37 ( 591) 


AMsLLE 


6 


-21.611 47 ( 562) 


1 X 10^ 


120 


+0.003 87 ( 572) 


AMsLLF 


12 


-75.778 67 ( 855) 


1 X 10^ 


431 


+0.014 28 ( 921) 


AMsLLG 


12 


-35.074 71 ( 683) 


1 X 10^ 


120 


-0.002 68 ( 697) 


AMsLLH 


6 


54.013 78 ( 619) 


1 X 10^ 


262 


+0.011 90 ( 664) 




12 


112.749 26 (1037) 


1 X 10^ 


512 


+0.007 52 (1071) 



dominated by Version A since Version B still does not have large statistics. Only subgroup 
IV(d) has been evaluated in both versions with comparable statistical weights 23| . Our best 
results for the gauge-invariant subgroups of group IV can be summarized as 



a 



a 



a 



a 



(8) 

IV{a) 
(8) 

IV{b) 
(8) 

IV{c) 
(8) 

IV{d) 



123.785 51 ( 44), 

- 0.417 04 (375), 
2.907 22 (444), 

- 4.432 43 ( 58), 



(58) 
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where 07^(6)' ^iv(c)i ^^"^ ^iv{d) statistical combinations of Version A and Version B. 

Summing up these terms one find that the contribution from all 180 diagrams of group 
IV is given by 

afi = 121.843 1 (59). (59) 
Finally, combining ipHjl with (jl^jl . (jTfj) and (jSHI), one obtains the value given in (fTH|) . 



VII. EVALUATION OF Af\m^/me,mf,/mr) 

There is a small contribution to from the three-mass term {rrif^/me, rrifj^/mr) which 
arises from 102 diagrams containing at least two closed fermion loops, of v-p and/or /-/type. 
The contribution of 30 diagrams analogous to those of Fig. ^ and Fig. |21 is 

Af^{m^/me,m^/m^) = 0.007 630 ( 1). (60) 

The contribution of 36 diagrams related to those of Fig. E]is 

Af^j{m^/me,m^/mr) =-0.053 818 (37). (61) 

The contribution of 36 diagrams analogous to those of Fig. is 

Af^y{m^/m,,m^/m^) = 0.083 782 (75). (62) 

Summation of these results leads to the value given in (fT7|) . The value 0.079 (3) quoted 



m 



24| is in rough agreement with (jU^ . In |24] it was assumed that the only nontrivial 



contribution to the eighth-order term arises from a muon vertex that contains an electron 
light-by-light scattering subdiagram and a tau vacuum-polarization loop and another in 
which the roles of electron and tau are interchanged. (See Fig. IHlwith (Zi, I2) = (e, r), (r, e).) 
It did not include the contributions ()6()j) and ()61|) . Our new calculation shows that this 
assumption was not justified. This is presumably because the mechanism that makes 
enhanced (see discussion in Appendix lA 2|) does not work if the momenta of photons ex- 
changed between muon and electron are not very small. 

Another term of order is 742^^(m^/mT-) which is calculable from 469 Feynman diagrams. 
However, its contribution to is of the order (m^/m,-)^ ln(mT-/m^)A2^''(l) ~ 0.005 so that 
it may be safely ignored without actual calculation. 
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VIII. DISCUSSION 



The size of integrals belonging to Groups I and II is rather small. Thus they have been 
evaluated using large number of sampling points, achieving precision of 5 or more digits. 
Furthermore, most of these integrals have been evaluated in alternative ways, either analytic 
or semi-analytic. The agreement between numerical and (semi) analytic calculations is so 
precise that it leaves no room for questioning the results. 

A similar comment applies to the integrals of Group III and Group IV(a), which are 
obtained by insertion of a vacuum-polarization loop 112 in the corresponding sixth-order 
diagrams, which have been fully tested against the analytic integration results. 

For the integrals of Group IV(d) formulation in Version B enabled us to discover an error 
in the Version A. After correcting the error, we now have two independent calculations 
which give the same results. For the remaining diagrams, of Groups IV(b) and IV(c), their 
structure had been tested extensively taking advantage of the fact that they have in general 
vertex and/or self-energy subtraction terms, which can be generated in two ways: One from 
a well-defined reduction procedure of the original unrenormalized integral, and another by 
construction of the renormalization terms from scratch, both analytically. The agreement of 
these two, proof of which often requires nontrivial analytic work, give a strong confirmation of 
their structure and of the master program from which all eighth-order integrals of individual 
diagrams is derived. See Appendix lA 21 for details. 

In order to obtain a further and definitive check, however, we have constructed IV(b) 
and IV(c) in Version B, too. As a consequence, we have integrands of Groups IV(b) and 
IV(c) in two versions. Extensive numerical work has shown that they are consistent with 
each other within the error bars of computation. This completes a comprehensive check of 
all diagrams contributing to Af\m^/me) by more than one independent methods. 

It is important to note that our classification of diagrams into groups (and subgroups) 
ensures that each subgroup is a gauge-invariant set. In fact, individual integrals are mostly 
not gauge-invariant and also infrared-divergent even when the UV divergence is renormalized 
away. On the other hand, each gauge-invariant set is well-defined and relatively small due to 
strong cancelation among its constituents. This is dramatically demonstrated by aiv{b) and 
ajv(c) in (EH), which are of order 1, whereas their constituents are two orders of magnitude 
larger as is seen from Table IXll 
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Empirically it is known that each mass-independent minimal gauge-invariant sets con- 
tributing to g-2 (namely diagrams containing no v~p loop or Z — / loop) is of order one apart 
from a power of a/vr. When v —p and/or / — / loops are inserted in such diagrams, they may 
acquire enhancements due to ln(m^/me) factor, which is a consequence of charge renormal- 
ization in case of v — p loop insertion and mass singularity present in the limit mg ^ in 
the case of / — / loop. The size of mass-dependent terms contributing to A2 {m^/rrie) may 
be understood semi-quantitatively from this observation. 

Let us now apply this argument to Mg2'2V2 Table |3 which is gauge-invariant and yet 
very large. Its size is inherited from the large sixth-order term Mqll. The extraordinary 



term with a large coefficient (6.38(8)) as 
It was noted then (unpublished) that 



size of Mqll is due to the presence of ln(m^/me~ 
was initially discovered by numerical integration 

it is numerically close to 27r^/3, later verified analytically j59|, enabling us to write 

27r2 

Mqll = — ln(m^/me) H . (63) 

Since Mqll is UV-finite, the term Inm^ comes from the scale set by the largest physical 
mass of the system, m^. The Inmg term arises from the integration of the momentum k of 
the I — I loop A4 over the domain Di {rrie < \k\ < nifj,, \pi\ < rrie), where pi, {i = 1,2,3) 
are the momenta of photons exchanged between the electron and the muon. Other domains 
such as D3(any k, \pi\ > rrie) does not contribute to Inm-e. 

What makes Mqll really large, however, is the presence of the coefficient vr^ in (jU^ . A 
physical interpretation for this fact was given by Elkhovskii jg^l who pointed out that, in 
the sub-domain D2 {riie < \k\ < m^, \pi\ « rrie, or more precisely \pi\ < anie, a ~ 1/137) 
the muon is nearly at rest and the electron can be treated as a non-relativistic particle in the 
field of the muon. One of the photons is responsible for the hyperfine spin-spin interaction, 
and the other two act essentially like a static Coulomb potential. It is the integration over 
the Coulomb photon momenta that gives a factor in each, contributing a factor 7r^(~ 10) in 

Actually it is not easy to maintain the nonrelativistic behavior of the electron throughout 
the domain Di outside of D2. This will result in the erosion of the enhancement factor vr^ 
in the domain Di — D2, although the ln(m^/me) behavior is still maintained. Together 
with non-logarithmic contribution from other parts of momentum space, the net effect is to 
reduce the contribution of the leading term of which is about 35, to about 21. This 
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reduction may be expressed crudely by choosing the fudge factor ^ to be about 0.12 in 

M^2r^) = ^ln(em,/m.). (64) 
Next consider the effect of the renormahzed photon propagator: 

D^^\q) = -^i^d^{q'/mla) + ■■■ , (65) 



where, to order a, 

dR{q^/ml,a) = 1 + - 

TT 



(66) 



When Dji is inserted in — 2 diagrams, the scale for the momentum q is set by the muon 
mass. This means that the leading term of the integral containing a vacuum-polarization 
loop 112 such as Mg2'2^p2 may be written as 

^6LL,P2 — "iKM^j^l^"^^ + terms linear in ln(m^/me), (67) 

where the factor 3 is the number of photon lines in which a vacuum-polarization loop can 
be inserted, and 

2 5 
K = - ln(m^/me) - g - 3> (68) 



provided that q^ is replaced by in ()66|) . The combination of these factors is responsible 
for the large size of the leading term of (|67p: 

3 X ((2/3) ln(m^/me) - 5/9) x 20 ~ 180, (69) 

which is 1.5 times larger than the calculated value of AfgLL P2 listed in Table 1X1 pretty close 
for a crude approximation. If the argument given below Eq. ()64|) is applicable to the photon 
momenta, we would obtain ~ 1.7 and (|69|) would become ~ 100. Both estimates would 
be acceptable as rough measure. 

It is important to note that Mg^'2''p2 is not only very large but also its value is known very 
precisely because it is obtained from the exactly known M^n by a well-understood vacuum- 
polarization insertion procedure. This means that the value of the term A2 {m^/nie) is 
determined primarily by ^/^(a) while its uncertainty comes mostly from 0/1/(6), a/v(c), and 
aiv{d)- This is why the value of A2\m^/me) did not change much even if a/y(d) suffered 
from a program error. 

These arguments may also be applied in identifying the leading terms of the tenth-order 
contribution A^2^\m^/me) Q|- 
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APPENDIX A: ELIMINATION OF ALGEBRAIC ERROR OF INTEGRALS 

Most eighth-order integrals, including those of Group IV, are huge. A systematic ap- 
proach is required to make sure that they are free from algebraic error and have forms 
suitable for numerical integration. To achieve this we adopted the following procedure: 

(a) Carry out momentum integration of Feynman diagrams and convert them into inte- 
grals over Feynman parameters using algebraic manipulation program such as FORM IgI]. 
This step is fully analytic. Conversion of all integrals of a gauge-invariant set is performed 
using a common template by permuting tensor indices of photon propagators. 

(b) Integrals thus obtained are divergent in general. Since computers are not capable 
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of handling divergence directly, both ultraviolet (UV) and infrared (IR) singularities must 
be removed from the integrand before integration is performed. In the subtractive on-shell 

n 

renormalization 16 2| of the n-th order diagram the renormalization term involving a m- 
th order vertex renormalization constant is of the form ~LmMn-mi where Mn-m is the 
g-2 term of order n — m. The subtraction procedure described in textbooks is not suitable 
for numerical integration, however, since it does not make the integrand of M„ — LmMn-m 
finite throughout the domain of integration, as long just a numerical constant. 

The first step to achieve a point- wise cancelation is to express as a parametric integral 
and combine it with the parametric integral of Mn-m using a generalization of Feynman's 
formula 



The domain of the combined integral ^ Mn-m may then be chosen to be identical with 
that of M„. Unfortunately, the integral is found to be intractable if we want to treat as 
a whole. However, if it is split as 



where the UV-divergent part Lm is identified by the highest power of f/, the IR divergent 
part Lm^'' by the highest power of V , Lm^'' Mn-m is found to have a term-by-term 
correspondence with UV divergent terms of the original (mother) integral M„, and UV 
divergences of M„ and Lm^^ ^ Mn-m cancel each other before (not after) integration is 
performed. (See ()A6|1 for the definitions of U and V.) IR divergence can be handled in a 
similar way. This point-wise subtraction is crucial for the success of our renormalization 
program on a computer. 

(c) In practice it is easier to start from the UV-divergent terms of the mother integral 
M„, which can be readily identified by a power counting at the singularity, and construct 
the subtraction term Lm^'' ^ Mn-m taking advantage of the term-by-term correspondence 
described above. This procedure, formalized as K-operation (see f)A32|) and ()A33|) ) can be 
applied to all orders. Further advantage of the K-operation is that it can be readily imple- 
mented in the FORM program that generates the integrand of M„. It is easy to confirm that 
the UV singularities of Lm'^^ Mn-m and M„ cancel out exactly by numerically evaluating 
the mother and daughter integrands at a sequence of points converging to the singular point. 




(Al) 



r _ r (UV) , r (IR) , r (finite) 
J^m — J^m ^ ^ ^m ) 



(A2) 
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(d) By construction the daughter integral factorizes analytically into a product of the 
divergent part of a renormahzation constant and a magnetic moment of lower order. Re- 
maining parts of renormahzation constants, such summed over all diagrams 
afterwards. The result is a convergent integral of lower order, which is easy to evaluate 
numerically. In other words, our renormahzation proceeds in two steps. But, of course, it is 
identical with the standard on-shell renormahzation. 

This procedure is designed to enable us to obtain extensive cross-checking at every step. 
To make this paper as self-contained as possible, let us describe them in some detail, although 
they can all be found in previous papers 0|. 



1. Construction of Feynman-parametric integral 

Let G be a 2n-th order proper lepton vertex of QED, which describes the scattering of 
an incoming lepton of momentum ^ — (^/2 by an external magnetic field into an outgoing 
lepton of momentum ^ + <^/2, where both leptons are on the mass shell. G consists of 2n 
lepton propagators and n photon propagators of the form (in Feynman gauge) 

+ ii + mi 
[ki + qiY - mj ' {ki + - mf ' 

besides factors describing the interaction, spinors, etc. Here ki is a linear combination of the 
loop momenta flowing through the line i. is a linear combination of external momenta, 
mj is the mass associated with the line which is temporarily distinguished from each other 
for technical reason. All these factors are combined to form a proper vertex part, which is 
integrated over n loop momenta. 

The first step of momentum integration is to replace + in the numerator of ()A3|) by 



an operator 



"OO o 

D^^^l dm\^ (A4) 



for each lepton line i. Since does not depend on ki explicitly, the numerator (turned 
into an operator now) can be pulled in front of momentum integration. The integrand then 
becomes just a product of denominators, which can be combined into one big function with 
the help of Feynman parameters zi, Z2, • • • , ^Af {N = 3n), assigned to respective propagators. 
Momentum integration can now be carried out exactly. can then be pulled back inside 
z integration. Omitting the factor (o/tt)" for simplicity the result can be expressed in the 
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form 



where 



r?») = ('^)"(„-i). /f.I^, (A5) 



N N N 



{dz)G = (^(l - ^i) n ^ ^^^^"i ~ ' ^'i' 



i=l i=l i=l 

N 



i=i ^ 
paud B., are homogeneous fonns of degree „ and „ - 1 in ...... . respectively. (See 

|57| for explicit definitions.) 
The operator has the form 

(X 7"n^i + ^i)7"'(^2 + m2)...7.---7"<'"-''(^(2n) +m(2„))7"f^"'. (AT) 

If G has closed lepton loops F;^ contains some trace operations, too. The action of F*^ on 
in (HSl) prod uces terms of the form 



V 



where the subscript k of stands for the number of contractions. By contraction we mean 
picking a pair Ipi + uii, pj + rrij from F*^, making the substitution 

(f. + m,, |)^. + m,)^(7^ 7^), (A9) 

multiplying the result with —^Bij, and summing them over all distinct pairs. Uncontracted 
Di are replaced by Q[. For k > 1, Fj^ includes an overall factor {n — 1)^^ ■ ■ ■ (n — k)~^ . 

In our problem it is convenient to use, instead of the vector Q/^ itself, a scalar function 
extracted from Q-^. Suppose — (external lepton momentum) enters the graph G at 
a point A, enters at a point G (which is the magnetic field vertex), and + ^7^/2 leaves 
at a point B. Then we may write 

Q;'^ = A^''\p>^ - + Af''\p^ + g72). (AlO) 

After a little manipulation we find, for example, 

1 ^ 

^AB) ^ ^AC) ^ ^CB) = --J2 riA^^B,, - P = PiAB). (All) 

i=i 
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associated with the path P = P{AB) of the external momentum p, which is any self-non- 
intersecting path starting at A and ending at B, and rjjp = (1, —1, 0) according to whether 
the line j lies (along, against, outside of) the path P. aI^'^\ etc., are called "scalar 

currents" since they satisfy an analog of KirchhofF's laws for electric currents when the 
diagram G is regarded as an electric network in which Feynman parameter Zi plays the role 
of resistance |6^]. Bij satisfies Kirchhoff's laws, too. U and Bij depend only on the topology 
of the graph G and not on whether the line is fermion or boson. They can be constructed 
easily from their definitions or by recursive relations starting from the one-loop case. We 
have written MAPLE and FORM programs to compute them algebraically for an arbitrary 
diagram. Once U and B^j are known, Ai = A^^^^ can be constructed by ()A11|) . For further 
details see js^]- 

The magnetic moment projection of rl?"'* of the muon in Version B is given by 

M^'"^^ = limTr[P"(p, q)Tl^% (A12) 

5=0 

where {p + = (p — if^Y = ^fi, and 

771 1 1 

P%P, q) = - 2^ + ^Mi^i - ii^p' - 3gV)(^ + 2^ + ^,)- (ai3) 

In the limit g = the term can be dropped in the denominator V of ()A5|) . Then V 
becomes a function of p^ only and can be simplified to 

V= z,m1-G, (A14) 

all leptons 



where G is defined by 

2"^ \dp 

G can be reduced further to the form 



G--\fi'Z\ ■ (A15) 



G= z.Aitnl, (A16) 

muon only 

by letting the external momentum p flows through consecutive muon lines only. This form is 
independent of how virtual photons are attached to muon lines. The information on photon 
attachment is contained in A^. This provides a significant simplification in programming. 
Eq. ()A12|) is the starting point of Version B. 
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In Version A the g — 2 term is projected out from the RHS of ()14p . In terms of Feynman 
parameters zi,Z2, . . . zn, {N = 3n — 1), introduced in the self-energy-like diagram G, the 
2n-th order magnetic moment can be written as 

M<-> ^ - 1). / (,.)(|±^j^ , iN . Z)^) , (A17) 

E, C, N, and Z are pieces of the magnetic projection defined by 

N = ^Tr[PM2GF)], E = \Tr[P,^E,], 

C = ^Tr[PrC,.], Z = irr[PrV]- (A18) 

The factors Pi and P2'^ , derived from the magnetic projector P*^ of ()A13|) by averaging over 
the directions of subject to the constraint q^p^ = 0, are of the form 



= 1 (1- + 1^ (g>^- - + —Y - ) • (A19) 



The operator F is the numerator factor of the self-energy-like diagram G: 

F = 7°H^i + mi)7"^(^2 + ma) . . . 7"1|)(2n-i) + m(2„_i))7"^ (A20) 

Here z and j refer to the internal photon lines arriving at the {2n — l) — th and 2r?,-th vertices 
along the lepton line (which depend on the diagram G), E'^ is defined by 

(9F 

E- = — = Yl ^^Fr, (A21) 

all leptons 

and FJ' is obtained from F by the substitution in the i-th line: 

ip, + rrii) =^ Y- (A22) 

Zfj_i, is defined by 

2n-l 

Z'^^ = ^ ;ZjZf , (A23) 
i=i 

where Z^'^ is obtained from F by the substitution 

{p, + m,) ^ \[YY{Pj + m,) - ip, + m,)Yr]. (A24) 
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C^"^ is defined by 



where 



C^ = J2C^JF^;^, (A25) 



i<j 



k=l l=k+l 

and F'fY is obtained from F by tlie substitution 



^ 2n-2 2n-l 

= 77^ E E ^''^liB'^kB-i - B'uB'^,), (A26) 



(^fc + mfe, pi + mi)^(Y, Y) (A27) 



in tlie k-th and /-tli lepton lines. See ()A6|) for tlie definition of B^|^. 

Note tliat the potentially troublesome in the denominator of ()A13|) is absent in the 
formula ()A18|) so that the limit q = can be taken without complication. As a consequence 
flA18|l depends only on one external momentum p, and the only scalar current needed are 
Ai of the muon lines associated with p. When Ai are expressed in terms of -Bjj's, they have 
the same form for all diagrams irrespective of how virtual photons are attached to the muon 
line. This provides a useful simplification in programming. 

We can now construct the integrand as follows: 

(I) Express the integrand as a function of symbols Bij, Ai, U, V for Version B and 
additional Cij for Version A. This can be achieved by an algebraic program (such as FORM 
01) by which momentum integration is carried out exactly. All integrals are generated from 
a small number of templates by permutation of tensor indices of photon propagators. 

(II) Quantities Ai, Bij, Cij, etc., introduced in (I) are just symbols. The next step is to 
turn them into explicit functions of Feynman parameters. This is facilitated by a common 
template which generates Bij and U for all diagrams sharing the same topological structure. 
Scalar currents Ai (and Bij) must satisfy up to eight junction laws and four loop laws for each 
diagram. This provides very strong constraints on scalar currents and sets up a powerful 
defense against trivial programming error. 

2. Construction of subtraction terms 

Following the discussion (c) at the beginning of this Appendix let us now discuss more 
explicitly how to construct UV-divergence subtraction terms starting from the mother inte- 
gral. 
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Suppose Mg is the magnetic moment contribution of a proper vertex part of 2n-th order 
defined by ()A5|) and ()A12|) . Carrying out the D-operations in Fj,, one finds 



, -^2 \ ('A281 

where mc is the maximum number of contractions of D operators, which is equal to n for 
a vertex part. In ()A28|) the Feynman cut-off of photon propagators is not shown exphcitly 
for simphcity. 

Suppose Mg has a UV divergence when all loop momenta of a subdiagram S consisting of 
Ns lines and ns closed loops go to infinity. In the parametric formulation, this corresponds 
to the vanishing of the denominator U when all Zi & S vanish simultaneously. (Note that Zi 
is a sort of conjugate to kf in a Laplace transform.) 

To find how a UV divergence arises from S, consider the part of the integration domain 
where Zi{E S) satisfy J2ies — ^- limit e — > 0, one finds 

V = t/ = C(e"^), 

B,, = 0(e"«-i) if i,jeS, 

= 0(e"«) otherwise. (A29) 

Let m5 be the maximum number of contractions of D operators within 5*. Simple power 
counting shows that the [m + l)-st term of Mg, whose numerator has at most ms factors 
of Bij, i,j G S, is divergent in the limit e ^ if and only if 

Ns — 2ns ^ min[m, ms], (A30) 

where min[m, m^] means the lesser of m and ms- 

If S" is a vertex part, we have Ns = 3ns and ms = ns- If 5* is an lepton self-energy 
part, we have Ns = 3ns — 1 and ms = ns — 1. In both cases, ()A30jl is satisfied only for 
min[m, ms] = ms, namely ms < m. This means that the UV divergence from S is restricted 
to the terms with ms contractions within 5* in the last mc — ms + 1 terms of ()A28|) . 

Let us now introduce a K5 operation, which extracts the UV-divergent part (due to the 
subdiagram S) of the Feynman integral 



Mg 
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= J{dz)GJG (A31) 



in an analytically factorizable manner. This is achieved by the following steps: 

(a) In the limit ()A29|) keep only terms with the lowest power of e in U, Bij, Ai, V. [U 
then factorizes into a product of Us and Ur, where Us, Ur are U- functions defined on S", R, 
respectively, and R = G/S is obtained from G by shrinking 5 to a point. Similarly for Bij. 
V is reduced to Vr. Factorization of Ai is more subtle since it has U in its denominator. 
The recipe here is to keep only those terms of Ai whose numerator is a linear combination 
of Bij with i,j G S. This is necessary for analytic factorization to work.] 

(b) Replace Vr obtained in (a) by Vr + Vs, where Vs is aV function defined on S. [Since 
Vs = 0{e) whereas Vr = this does not affect the leading singularity of the integrand 
in the e — limit. 

(c) Rewrite Jq in terms of the redefined parametric functions, drop all terms except those 
with ms contractions within S, and call the result K.sMg- 

Since we deal in practice with logarithmic divergence only, the steps (a), (b) and (c) are 
sufficient to ensure that (1 — K.s)Mg is convergent for e ^ 0. The inclusion of Vs in (b) 
serves two purposes. One is to avoid spurious singularity which Vr alone might develop in 
other parts of the integration domain, and the other is to enables us to decompose K^Mgi 
(assuming 5* is a vertex part) into a product of lower order factors analytically as 

KsMg = LsMg/s, (A32) 

where Ls is the UV-divergent part of the renormalization constant Ls- 

If 5* is a lepton self-energy part inserted between consecutive lepton lines i and j of G, 
we obtain a slightly more complicated result 

KsMg = 6msMG/s(^*) + BsMG/s(i') (A33) 

Here i3s and 5rhs are UV-divergent parts of renormalization constants Bs and 5ms. Since 
Ls,Mg/s, etc. are quantities of lower-orders, they are already known or can be easily 
constructed from scratch. 

Note that Eqs. ()A32|) and ()A33|1 are quite nontrivial since the LHS are defined over the 
entire n-dimensional surface while the RHS are products of integrals over lower- dimensional 
spaces. Identification of the second term on the RHS of ()A33|) requires further work involving 



an integration by part. (See [57[ for details.) 

An IR divergence, which is caused by vanishing virtual photon momentum, arises from 
the part of integration domain R where Zi's assigned to the photon takes the largest possible 
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value under the constraint Yl^i — 1- other Zj's are pushed to zero in the IR hmit. 
Furthermore, the IR singularity, in order that it actually becomes divergent, must be en- 
hanced by vanishing denominators of muon propagators adjacent to the infrared photon. 
In parametric language this corresponds to the vanishing of V a.s V ^ 6"^ for 5 ^ in the 
integration domain characterized by 



Starting from this one can obtain a power counting rule which enables us to identify IR 
divergent terms in a way analogous to that of UV divergence. The criteria to be satisfied 
by the IR subtraction operator 1^ are 

(i) it subtracts the IR divergent part of the mother integrand completely, 

(ii) it factorizes analytically into a product of IR-divergent part of renormalization con- 
stant and magnetic moment of lower orders. 

The difference with K5 operation is that we must now look for largest negative powers 
of V instead of U in ()A28jl . See for details. 

3. Diagrams containing a second-order lepton self-energy part 

When an integrand containing a second-order electron self-energy part 5* is expressed as 
a function of scalar currents, all of its terms contribute to the UV divergence. This means 
that the integrand of the self-energy subtraction term, when expressed in terms of its own 
scalar currents, has a form identical with that of mother integrand. Their difference comes 
solely from different structures of scalar currents for mother and daughter integrals. This 
provides the simplest example of ()A33|) . 

To demonstrate it explicitly let us go back to the vertex ()A5|) and rewrite F,y of ()A7|1 to 
exhibit the self-energy part S explicitly: 



0{6) if i is a muon line in R 



Zi = \ 0(1) if i is a photon line in R 
0{6^) if i is in G/R. 



F 



X (i+i) + m(i+i)) . . . 7"(2"-i) (^(2,) + m(2„))7' 



(A34) 
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where 7^(^(j) + m(^i~^)'-^p comes from the second-order lepton self-energy part S. {jy is not 
shown exphcitly.) Carrying out the contraction of 7'^ and 7/3 this factor can be written as 
2m(i) — 2(_^(j) — m(j)). This leads naturally to the decomposition 

F. = F«+F(2), (A35) 

where F^,^'' and F^^^ are obtained by replacing 7^(_p(j) +m(j))7^ of ()A34|) 2m(j) and — 2(_p(j) — 
m(j)), respectively. 

As is easily seen the fI^^ part of the integral K5MG factorizes exactly into a product of 
the self- mass 5m2 and the term MG/s{i*)i which is obtained from Mq by shrinking the lepton 
self-energy diagram S* to a point: 

5m2MG/s{i*). (A36) 
In the same limit the F^. part of the integral K5MG factorizes exactly into a product of 



B2 (the UV-divergent part of B2 (see 0])) and the term M^^gf^^i^, which is obtained by 
inserting the factor — rn(^j;j) in Mc/sii*)'- 

B2M^/s(,'y (A37) 

Note that in the K^Mg operation leading to ()A37|) the contraction of with other D's 
in G are dropped. 

Finally, using an identity 

which is a particular case of Eq. (A. 7) of Nakanishi's Appendix 65|, one can transform 
^G/s{i) amplitude MG/[s,i+i], which is obtained by removing the self-energy part S 

and the adjacent lepton line z + 1 from G. (This is identical with a vertex of lower order 
obtained directly by Feynman-Dyson rules.) 

4. An Illustration: Two-Step Renormalization of Fourth-Order Magnetic Moment 

The formulation described above is illustrated here by applying it to the fourth-order 
magnetic moment M4 in Version A, which consists of two parts M^a and M^b. Eq. ()A17jl 
applied to the diagram G = [1, 2, 3, 4, 5] of Fig. IHl^d) leads to 

M.a = -^j (dz) + + J , (A39) 
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FIG. 8: Assignment of Feynman parameters zi . . . is shown for (a) second-order self-energy 
diagram, (b) second-order vertex diagram, (c) second-order self-energy diagram in which a 2-point 
vertex is inserted, (d) fourth-order diagram M4a, (e) fourth-order diagrams M4f,. Horizontal lines 
(except in (b)) are lepton lines in the presence of the magnetic field. Each of diagrams (d) and (e) 
represents the sum of three fourth-order vertex diagrams. 

where , for simplicity, Feynman cut-off of photon lines is not shown explicitly [stI . Is? ]. 
Numerator functions are expressed in terms of scalar currents Ai and Bij: 

Eo = 8{2A,A2A; - A,A2 - A.A^ - A^A^), Co = -24z4Z,/U, 
No = G{Eo-8{2A2-l)), 

Zo = 8zi{-Ai + A2 + A3 + A1A2 + A1A3- A2A3) 

+ 82:2(1 - A1A2 + A1A3 - A2A3 + 2A1A2A3) 

+ SzsiAi + A2-A3- A1A2 + + 
TVi = 8G'[fii2(2-A3) + 5i3(2-4A2) + 523(2- A)], 

Zi = -8^i[5i2(l- A3) +5x3 + 523^1] 

+ 8Z2[B^2{1 - A3) - 4S13A2 + ^23(1 - Ai)] 

- 8^3[5i2A3 + 5i3 + ^23(1 - Ai)], (A40) 
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in addition to zi,Z2,zs, and G, and 

Bn = 2:235) = z^^, Bi^ = —Z2, -B23 = ^u? 

B22 = 2:1345, i?33 = 2:124, U = Z2B12 + Zi4-Bii, 

Ai = 1- {ziBu + 2:2^2^ + zA)/U, i = 1, 2, 3, 
G = z^A^ + Z2A2 + 2:3^3, V = 2:123 - G, 
{dz) = (5(1 - 2;i2345)c?2:iC?2:2C?2:3C?2:4rf2:5, 2:235 = 2:2 + 2:3 + 2:5 



(A41) 



The integral for M^b of Fig. IHl^e) has the same form as ()A39j) but has different definitions 



of functions 
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= 8Ai[4(A2 - Ai) - A1A2I Go = -8A2, 

No = -8G[4(1 -A, + Al) + ^2(1 - 4Ai + Aj)]. 

Zo = 8z,3[AA,-A2{l + Al)]+8z2A2il + Al), 

iVi = 8G[8(5ii-fii2) + 3AiEi2], 

Zi = 24(^13 - Z2)AiBu, 



(A42) 



where 



Bu — ^24, B12 — Z4, , S22 — ^1345? U — Z135-B11 + 2:2-812, 

Al = Z55ii/f/, A2 = z.Bu/U, G = z,sA, + 22^2, V = zi23 " G. (A43) 
The standard on-shell renormahzed amphtudes a^a and 04;, are defined by 

aia = Mia - 2L2M2 (A44) 

and 

046 = - 5m2M2* - B2M2. (A45) 

We carry out the renormalization in two steps, expressing all quantities as parametric inte- 
grals. For instance, the magnetic moment M2 is written in the form (putting = in 



where 



M2 = -\j {dz)^^^^, No + Zo = 4G(Ai - 1), (A46) 



Bii = 1, U = zi4,Bii, Al = Z4,/U, G = ziAi, V = zi + m\z/i — G, 

[dz^ = 5(1 — Ziijdzidzi, Zii = 2:1 + 2:4. (A47) 
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Following the general discussion the parametric integrals of B2 and Sm2 are split into 
UV-divergent parts B2 and 6m2 and UV-finite part B2 and 67712: 



B2 = B2 + B2, 5m2 = 5m2 + Sm2, (A48) 



where 



B2 = - J (dz) J^^ zdmi-^, Eo = -2A,, 
B2 = \j(dz)^, iVo = 2G(4-2AO, 

5m2 = \j {dz) j^^ z^dml-^, Fo = 2(2 - Ai), 5m2 = 0. (A49) 

Ai, U, V, (dz) are the same as in ()A40|) and A and A are UV and IR cut-offs for the virtual 
photon mass m^. The UV cutoff is removed from B2 since it is not UV-divergent. 

Similarly, we split the vertex renormalization constant L2 for Fig. Efb) as L2 + L2 , where 

1 f f^^ F 
L2 = --J (dz) J ^ Zidmljj^, Fi = -2Bii, 



^^{dz)jj^, Fo = -2(1 -AA, + Ai), (A50) 
where 

Bu = 1, U = zi2aBii, Ai = z^/U, G = z^Ai, V = zu + m\zA - G, 

(dz) = (5(1 — Zi24:)dzidz2dz4. (A51) 

Finally we need M2* (magnetic moment contribution of the diagrams in Fig. IH^c)): 

N* + Z* = -8GAi{Ai - 1), A^* = IQGBn, E* = -8Al, (A52) 
where functions are the same as in ()A51|) . 

(a) Two-step renormalization of M^^a- 
The first step is to rewrite ()A44|) as 

a4a = AM4a + ((K12 + K23)M4a - 2L2M2). (A53) 

K12 is a UV-divergence extraction operator for the vertex subdiagram S = [1,2,4]. The 
integral 

AM4a = (1 - Ki2 - K23)M4a (A54) 
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is UV-finite by construction. K.uMia can be written as 

1 r f^^ n' + z' 

K12M4. = -j (dz) J^^ zdml ^,3"^,/ , (A55) 
where the Feynman cut-off of photon mass is shown exphcitly and 

-B12 = ^35) U = Zi2iBi2, A^ = Zrj/z^^, A^=Z^lzx2^ 
G' = Z^A'^, V' = Zi23 + mlz4 - G' - ^12^1, 

N[ + Z[ = 8G'Bi2{1-A'^). 

K23M4a for the vertex diagram [2,3,5] can be constructed in a similar way. 

To show that ()A55j) can be factorized and reduced to the RHS of ()A32j) . let T = [3,5] 
be the reduced diagram obtained from G by shrinking 5* to a point. Let us define zs = 
Z124, zt = -235, and introduce Feynman parameters Xj and yj for the sets S and T in such a 
way that 

xi = zi/zs, X2 = Z2/ZS, X4 = ^4/2:5, Vs = zs/zt, 2/5 = z^/zt- 
Then, after few steps of manipulation, which amounts to the substitution 

Bu^Bl2{=l), A',^Al{=y,/y,,), G' ^ y^Al 

U'^UsUt, V' ^zsVs + ztVt, N[ + Z[ ^ 8y,Al{l - AJ) , 

with 

Af = X4/X124, Us = zs, Ut = Zt, 
Vs = X12 + m\xA - Gs, Vt = yz~ Gt, Gs = xuAf, Gt = 2/3^3 > 
we can rewrite ()A55|) as 

KM. = ^ / / m I (rfj) .,<;,„jg|il_0. (A56) 

where (dz) = 5(1 — zs — ZT)dzsdzT- Note that U' in ()A55|) factorized as UsUt, and they 
canceled out completely against zs and zt in the numerator. Since zs and zt now appear 
linearly in the denominator only, the z-integration can be carried out explicitly using Eq. 
()A1|) which reduces the integral ()A56|) to a product of L2 and M2: 

KuM^a = L2M2. (A57) 
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This is a particular case of ()A32|) . Similarly for K23M4a. Note that K-operation extracts 
only the UV-divergent part L2 of L2. This exact factorization of ()A57|) is generalizable to 
all orders, which is crucial for carrying the two-step renormalization scheme over to higher 
orders. 

From (|Xi4ll . (IX50|l . and (1X571) we obtain 

a4a = AM4, - 2L2M2, (A58) 
where both terms on the right-hand-side are UV-finite. (L2 is IR- divergent.) 

(h) Two-step renormalization of M^,. 
We begin by rewriting ()A45|) as 

= A'M4i, + (K2M4fe - 5m2M2. - B2M2). (A59) 

K2 isolates the UV divergence of M^i, arising from the self-energy subdiagram S = [2,4]. 
By construction the integral 

A'Mih = (1 - K2)M4f, (A60) 



is UV-finite. Following the rule given in Appendix IA21 'K2M^b is obtained by dropping 
terms in Zq and Zi of ()A42|) that contains Z2 explicitly {z2 and Z4 hidden within Ai and Bij 
must be kept to the leading order in the K2 limit.) Noting that A2 — > A'^A^ in the K2-limit 
we find from ()A42|) 

= 8A['[4{A', - 1) - A[A',], = -8A[A'„ 
iV; = -8G'[4(1 - a; + A[^) + A,A'^{1 - AA[ + A[% 

= 8z,sA[[A-A',il + A[% 
N[ = 8G'[8iB[,-B[,) + 3A[B[,], 

Z[ = 2AzisA:,B[^, (A61) 

where 

Ai=Z5/zi35, ^2 = ^4/^24, G =^13^4^, 

Including the regulator mass of photon 4 explicitly, we can express K2M4b as 



K 



'Mib = Yq j ('^^)^4 j dml 



K + C'o 2{N', + Z;) N[ + Z[ 
U'2Y'2 + u"^V'^ U'^V'^ 



(A62) 
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In the K2-limit U' and V' decompose as 

U' U'sU't, V' zsV's + ztVt, (A63) 

where U'^ = zs = Z24, = zt = -2135, and Vg and are V functions for the sub diagrams 
S = [2, 4] and T = [1, 3, 5] to be made exphcit in the following. 

Let us introduce Feynman parameters Xi and yi for the subdiagrams S and T as follows: 

X2 = Z2/zs, X4 = Z4/zs, yi = Zi/zT, ?/3 = ^3/^T, y-> = Z^/zT, 

and rewrite (jA62|l in terms of 

= X24ZS, B[2 = x^zs, A[ = 1/5/1/135, A'2 = X4/X24, Gs = Gt = yisK: 

U's = zsU's, U's = X2i, U't = ztU'^, U'^ = 1/135, 
V' = zsVg + ztV^ , Vs = X2 + mlx4-Gs, V^t = 2/i3 - ^'r. 



(dx) = S{1 - X24)dx2dx4, (dy) = 5{1 - yi35)dyidy3dy5. (A64) 
Then, dropping superscripts ' and " for simplicity, we can re-express ()A62|1 as 

K2M46 = ^ y i'^^) J {^y) j i.^^) j Xidm\ 

. zlr^z^-^ SAl{A{A2 - 1) - A1A2) - 8A1A2 

4"'4"' -16Gt(4(1 - Ai + Al) + ^1^2(1 - 4Ai + Al) - (4 - ^2(1 + Al))) 



+ 



Up^ {zsVs + zrVrf 



zf^z^-^ 8Gt{8{Bu - B12) + 3Ai5i2 + 35 



where 



(c?5) = (5(1 — 25 — z^jdzgdzT- 
Now the 5 integration can be carried out exactly using ()A1|) . The result is 



+ 



IQJ J Jx^ UIVsU^Vt 
-8G't(4(1 - a + Al) + ^1^2(1 - 4Ai + Af) - (4 - ^2(1 + Af))) 

UlVsU^VS 



, 8GTfiii(8(l-A2)+3AiA2 + 3A2). 
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where B12/B11 = A2 is used. Decomposing the numerator of each term into terms propor- 
tional to (2 — A2) and A2 following ()A35|) . we can rewrite ()A66|) as 



K2M4b 



+ 



1 

16 

X 



16 



X 



(dx) 

(dy) 
(dx) 

{dy)[ 



-16AI -16Gt{-Ai + Af) 



A2 



32,(3x1311 
U^Vt 



x^dml- 

A2 f^S^S 



A2 



8^2(2 - Ai) - 8A1 ^ -8Gt(1 -Ai-Aj + Al) 



+ 



U^Vt 
2G'rEn(-4 + 12Ai) 



(A67) 



Note that both terms are now just products of x-integral and y-integral. Comparing them 
with (fX49|l and |X52|l . it is easy to see that the first term is 5m2M2*, which corresponds to 
the first term of ()A33|) . The x-integral of the second term is proportional to -B2. Thus the 
remaining task is to identify the y-integral with M2 in ()A49|) . The first step of demonstration 
is to transform M2 by an integration-by-part. Dropping the suffix T for simplicity, we obtain 
WiAi(l-Ai) 



Mo 



(dyy 



{dy)yi 

where we have used 



y^A^{l-A,] 
A(l-Ai) , y^Bn 



y,A,{l-A,){l-Al] 



(A68) 



B 



11 



1, U = yuBn, A, 



y^B 



11 



dU 
dyi 



B 



dAi 



AiB 



11 



ii> 



U ' dyi dyi 
Noting that yiBu/U = 1 - Ai, (|X68|l can be reduced to 



U 



^-1-A^ 



Mo 



{dy)yi 



2Ai(l - A^){l - 2Ai) ^ yiA,{l - A,){1 - Aj) 



(A69) 



It is now easy to see that (jA69j) is proportional to the second y-integral of ()A67jl . One 
has only to note that {dy)yi = 6{1 — yi — yi)yidyidy4^ in ()A69jl is equivalent to {dy) = 
5(1 — Vi — ys — y5)dyidyzdy^ in ()A67jl since its integrand depends on the sum yi + y^ only. 
Eq. (|X68|l is a special case of ()A38|) . 

Making use of formulas for 6m2, B2, and M2* in ()A49|) and ()A52|) . we thus obtain 



K2M46 = 5m2M2* + B2M2, 



(A70) 
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which is a particular case of ()A33|) . 

Using (IX59|l and (IX60|) ciAb ill ()A45|1 can be rewritten as 

= A'M4b - B2M2. (ATI) 

Both terms on the right-hand-side are UV- finite although A' M41, is IR-divergent. 

(c) Separation of IR divergence in A' M41,. 

The sum a^a + ^Ab is UV- and IR-finite. However, for numerical evaluation it is necessary 
to remove IR-divergent terms from the integral A' M^h- This is where Step 2 comes in. A 



general procedure for isolating the IR divergent terms is given in [57[ . Power-counting shows 
that, of three vertex diagrams contributing to A' M^i,, IR divergence arises only from the 
vertex in which the magnetic field acts on the muon line denoted Z2 in Fig. ^e). More 
specifically, it arises from the domain 

^5 = 1 - 0{6), z,, Z3 = 0{6), Z2, z, = 0{6'), (A72) 

where 5 ~ A and A is the photon mass. In this domain one may define an IR-extracting 
operator It by 

U^UsUt, V^Vs + Vt, Vs = Z2{l-A2), Vt = z^il - A^) + z^ml 

Ai = Z5/zi35, A2 = Z4/Z2A, Us = Z24, Ut = Z135. (A73) 

This definition is actually identical with that of the K2-limit ()A62|) except that the 
substitutionA2 A'^A'2 (see above Eq. ()A61|) ) is replaced by A2 A^i^ A2). 
The separation of IR-divergent term of A' M41, may be written as 

A'M46 = AM4b + IrA'M^b, 

where 

AM46 = (1 - lT)A'M4fe. (A74) 

The definition ()A73|) . although it picks up the IR-singularity correctly, is actually not 
complete until an appropriate numerator function is chosen. We chose to define the IR 
separation operator It, T = [1, 3, 5], by 

lTA'M4, = ^|(rfz)^, (A75) 
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where 

Ft = -2(1 -4^1 + ^2), Fs = -4^2^2(1-^2). (A76) 

Ft is chosen to coincide with Fq of ()A50|) and Fs corresponds to Nq + Zq of ()A46|) . After 
taking steps analogous to ()A63|) and ()A64|) . this choice leads to an exact factorization of the 
integral ()A75j) as a product of known second-order integrals 

IrA'M^b = L2M2. (A77) 

Note that a particular choice of numerator is not crucial as far as it deals with the IR 
divergence correctly since what is subtracted must be put back in the end. 
Summing up all terms we obtain 

= AMia + AM4b - AB2M2 (A78) 

where AB2 = L2 + B2{= 3/4) and M2 = 1/2. Numerical evaluation of AM^a and AM4b 
gives 

AM^a = 0.218 342 (17), 
AMib = -0.187 501 (14), 

04 = -0.344 158 (22), (A79) 

which is an update of the old evaluation [57 1. It is in good agreement with the analytical 
value -0.344 166 ... . 

We described the fourth-order case in full detail because it will serve as a good prototype 
for higher-order cases. To begin with integrals such as lA39j) with integrands ()A40jl and ()A42|1 
are obtained by a simple algebraic program written in FORM. Subsequent manipulation of 
integrands proceeds in a well-organized manner. The important point is that all higher-order 
integrals can be handled in the same manner, the necessary extension being straightforward. 
This is the reason why we are able to treat the algebra of higher-order cases with complete 
confidence. 
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APPENDIX B: NON-STATISTICAL ERROR IN NUMERICAL EVALUATION 
OF FEYNMAN INTEGRAL 



Our integrand of Group IV is an algebraic function of more than 4,000 terms, each term 
being a product of up to 10 functions defined on a unit 10- dimensional cube: 

0<Xi<l, i = 1,2,--- ,10. (Bl) 

FORTRAN codes of some integrals are as large as 100 kilobytes. These integrals are identical 
with the corresponding integrals for the electron vertices, the only difference being the value 
of the parameter me/m^. However, the behavior of muon integrals are strongly influenced 
by the presence of a singularity located at a distance of order (me/m^)^ just outside of 
the integration domain ()B1|) . i. e., a unit cube. This makes numerical integration of some 
integrals more delicate or difficult compared with the electron case. This is the main (though 
not the only) source of the d-d problem in the muon g — 2 calculation. 

Numerical integration of these integrals is carried out using an adaptive-iterative Monte- 
Carlo integration routine VEGAS . It is the only effective method currently available to 
integrate such huge integrals. It is an adaptive-iterative integration routine based on random 
sampling of the integrand. In the i-th iteration, the integral is evaluated by sampling it at 
points chosen randomly according to a distribution pj_i (a step function defined by grids) 
constructed in the (z — l)-st iteration. This generates an approximate value of the integral, 
its uncertainty cij, and the new distribution function pi to be used in the next iteration. The 
distribution pi is constructed in such a way that the grids concentrate in the region where 
the integrand is large. The construction of pi in the {i — l)-st iteration involves a positive 
parameter f3 that controls the speed of "convergence" to a stable configuration. In most 
cases we chose P = 0.5. We may even be forced to choose /5 = (no change in p), which is 
necessary in some difficult cases. 

After several iterations Jj and (Xj are combined assuming that all iterations are statistically 
independent. The combined value and uncertainty are given by 

/ = (E(^^/^'))/(E(iM')), ^ = (E(v^n)^'/'- (B2) 

i i i 

For well-behaved integrals pi converges rapidly to a (practically) stable configuration. Once 
Pi is stabilized, the error generated by VEGAS is (nearly) statistical and proportional to 
A/"^^/^, where M is the total number of data samplings. 
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After the point-by-point renormalization is made the integrand has the form 



/ = /o + ■ ■ ■ + /, 



(B3) 



where /o is obtained directly from a Feynman diagram and /i,..., fr are terms needed to 
renormahze UV (and/or) IR divergences of /q. Terms /o, fr are all divergent on the 
surface of the unit cube ()B1|) . The sum / is mathematically well-defined and integrable. 

This does not guarantee, however, that / is well-behaved on a computer. This is because 
expressions for /o,..., fr on computer are only as accurate as the number of digits in use (64 
bits, 128 bits, etc.). In the part of the domain where /o, fr are singular, / loses most or all 
of significant digits and is affected severely by round-off errors. When this happens, Jj and 
(Tj become unreliable or even divergent. Note that this problem is an inevitable consequence 
of any computer calculation in which only a finite number of digits is available. We shall 
call it the digit-deficiency or d-d problem. In order to cope with the d-d problem before it 
upsets the integration, we have developed several strategies. 

a. Stretching. The integrand / defined in ()B3|) may still have integrable singularities on some 
boundary surfaces, which can be removed by an appropriate change of variables. However, 
it is difficult to find analytically correct mapping because of the complicated structure of the 
integrand. A simple way to remove or weaken the d-d problem is the "stretching" defined 
as follows: Suppose VEGAS finds after several iterations that the integrand samplings tend 
to concentrate in the vicinity of an (n — l)-dimensional surface, say xi = 0, perpendicular 
to the xi axis. Then the mapping 



where ai is a real number greater than 1, stretches out the domain near xi = and random 
samplings in the x[ variable give more attention to this region from the beginning of iteration. 
Also, the Jacobian aix"'^~^ of this mapping weakens the singularity. Similarly, the singularity 
at xi = 1 can be weakened by the stretching 



Stretching is a one-to-one mapping of a unit hypercube onto itself. It may be applied to all 
variables independently. An appropriate stretching speeds up convergence of p considerably. 
Note also that different stretchings lead to statistically independent samplings of an integral 




(B4) 




(B5) 
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which must give the same answer within error bars. This flexibihty is important in assessing 
the rehabihty of results of integration. Of course, stretching does not always work well since 
it disregards the actual (and hard-to-identify) analytic structure of the integrand. 

b. Splitting. Going from double precision (real*8) to quadruple precision (real* 16) (or 
even higher) arithmetic is the most effective way to control the d-d problem. One practical 
obstacle is that real*16 slows down computation by a factor 20 ~ 30. Thus we were not 
able to use real* 16 extensively until massively parallel computers became readily available. 

Actually, in many cases, real*16 is needed only in a small part of the integration domain. 
It is therefore useful to adopt the following strategy: Start the evaluation of a Feynman 
integral in real*8, which explores the integrand at high speed. If it identifies the region 
causing the d-d problem, split the integration domain into a small (rectangular) part in 
which the d-d problem occurs and the remainder. The difficult region is then evaluated in 
real* 16, while the rest continues in real*8. This strategy has been very successful and most 
integrals have been evaluated in this manner. 

Recently, a modified algorithm of VEGAS has been developed which makes this splitting 
local and automatic j^]. In this approach the integrand / is first evaluated at each point 
in real* 8. The result is tested by computing the ratio 

t = {U + \f-\)/\U + f-l (B6) 

where fj^{fJ) is the sum of positive (negative) terms of /o, ■■■fr- If ^ is larger than an empiri- 
cally selected number to? it signals a possible d-d problem. The integrand is then reevaluated 
in real*16 at the same spot. If the d-d problem is not severe, this method is very efficient 
and runs much faster than pure real*16. In difficult cases, however, a simple splitting may 
work faster since it does not require the overhead needed in computing ()B6|) . 

c. Freezing. Sometimes, it is very difficult to find a reasonable stretching that does not run 
into the d-d problem before it settles down to a (nearly) stable p. In such a case, one may 
freeze p by putting /? = few steps before the d-d problem becomes serious. The resulting 
p is not optimal so that it requires longer computing hours to achieve the desired statistical 
uncertainty. 

d. Chopping. If procedures a, 6, c fail to solve the d-d problem, one may restrict some 
integration axis (0, 1) to (5,1 — 5), where < 5 ^ 1, to exclude the danger zone. This 
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is referred to as chopping. The error introduced by chopping is of order S^^'^(\nS)'^, where 
the positive number a can not be fixed without knowing the analytic structure. In practice 
it is sufficient to find a crude value of a by carrying out integration for several values of 
6. When chopping is used, we must carry out full scale calculations for several 6 (which 
require extra computing time). Note also that integration becomes more and more difficult 
as 6 gets smaller. The difficulty in assessing the effect of chopping was the major source of 
non-statistical uncertainty in earlier calculations. 

Chopping can produce a crude approximate result vary rapidly. Thus it was used in early 
stage of our work to obtain estimates of a rough order of magnitude. However, it turned 
out to be not effective for obtaining more precise results. Thus it was abandoned entirely in 
the later phases of our work. 

Our final results were obtained using stretching, splitting, and freezing, or their com- 
binations. In most cases stretching and splitting are sufficient to solve the d-d problem. In 
some cases, however, even splitting was not sufficient. In the absence of higher precision 
arithmetic, the only effective way to control the d-d problem was freezing. In such a case 
it is still useful to divide the integration domain into several pieces and apply freezing in 
only one of them. Moreover, it is found to be useful to restrict the number of samplings per 
iteration to a relatively small number and use a very large number of iterations. This will 
enable us to accumulate large statistics while controlling the amount of wasted iterations to 
an acceptable level. 
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